// INITSPIN.C
// Functions for initialization of spins.
// © 2021 Peter J. Meyer
#include "iss.h"
#define ISING_UP_SPIN_PROBABILITY \
((num_spins*(1+initial_magnetization)-2*num_up_spins) \
/(2*(num_spins-num_up_spins-num_down_spins)))
#define Q_STATE_POTTS_UP_SPIN_PROBABILITY \
(((double)num_spins/q)*((q-1)*initial_magnetization+1)-num_up_spins) \
/( num_spins-num_up_spins-num_other_spins)
static void assign_initial_ising_spins(void);
static void assign_initial_q_state_potts_spins(void);
// On return from this function abs_magnetization
// and abs_internal_energy have been calculated.
/*---------------------------*/
void assign_initial_spins(void)
{
if ( model_is_ising )
assign_initial_ising_spins();
else // model is q-state Potts
assign_initial_q_state_potts_spins();
}
// Call with initial_magnetization = 1.0 for all spins up,
// -1.0 for all spins down, 0.0 for random up- and down-spins.
// This function does not assume any particular state of the spins.
// It also calculates the absolute magnetization.
/*----------------------------------------*/
static void assign_initial_ising_spins(void)
{
int i, j, k, l;
int num_up_spins=0;
int num_down_spins=0;
switch ( dimensionality )
{
case 2:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
if ( site_diluted ) // site_diluted set in SETPARAM.C.
if ( !SITE_IS_OCCUPIED_2(i,j) )
continue;
// Site is occupied.
if ( initial_magnetization == 1.0 )
{
sites2[i][j] &= BIT0_0; // Set spin to up.
num_up_spins++;
}
else if ( initial_magnetization == -1.0 )
{
sites2[i][j] |= BIT1_0; // Set spin to down.
num_down_spins++;
}
else
{
// Under some (rather unlikely) scenarios ISING_UP_SPIN_PROBABILITY
// could be < 0 or > 1, but this is not a problem.
if ( PRNG1() < ISING_UP_SPIN_PROBABILITY )
{
sites2[i][j] &= BIT0_0; // Set spin to up.
num_up_spins++;
}
else
{
sites2[i][j] |= BIT1_0; // Set spin to down.
num_down_spins++;
}
}
}
}
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;
// Site is occupied.
if ( initial_magnetization == 1.0 )
{
sites3[i][j][k] &= BIT0_0; // Set spin to up.
num_up_spins++;
}
else if ( initial_magnetization == -1.0 )
{
sites3[i][j][k] |= BIT1_0; // Set spin to down.
num_down_spins++;
}
else
{
if ( PRNG1() < ISING_UP_SPIN_PROBABILITY )
{
sites3[i][j][k] &= BIT0_0; // Set spin to up.
num_up_spins++;
}
else
{
sites3[i][j][k] |= BIT1_0; // Set spin to down.
num_down_spins++;
}
}
}
}
}
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;
// Site is occupied.
if ( initial_magnetization == 1.0 )
{
sites4[i][j][k][l] &= BIT0_0; // Set spin to up.
num_up_spins++;
}
else if ( initial_magnetization == -1.0 )
{
sites4[i][j][k][l] |= BIT1_0; // Set spin to down.
num_down_spins++;
}
else
{
if ( PRNG1() < ISING_UP_SPIN_PROBABILITY )
{
sites4[i][j][k][l] &= BIT0_0; // Set spin to up.
num_up_spins++;
}
else
{
sites4[i][j][k][l] |= BIT1_0; // Set spin to down.
num_down_spins++;
}
}
}
}
}
}
break;
}
if ( adjust_zero_initial_magnetization )
if ( num_up_spins == num_down_spins )
{
flip_initial_spin_to_up();
num_up_spins++;
num_down_spins--;
}
abs_magnetization = num_up_spins - num_down_spins;
if ( internal_energy_measured )
// Since the bonds have already been established we can now:
abs_internal_energy = get_abs_internal_energy();
}
// Call with initial_magnetization = 1.0 for all spins up,
// -1/(q-1) for all spins in another direction, 0.0 for random spins.
// This function does not assume any particular state of the spins.
// It returns the number of up spins minus the number of other spins.
/*------------------------------------------------*/
static void assign_initial_q_state_potts_spins(void)
{
int i, j, k, l;
int spin_value;
int num_up_spins=0, num_other_spins=0;
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;
// Site is occupied.
sites2[i][j] &= ~spin_value_mask; // Set spin to up.
if ( initial_magnetization == 1.0 )
// Spin is up.
num_up_spins++;
else if ( initial_magnetization == -1.0/(q-1))
{
// Set spin to another direction.
spin_value = 2 + (int)(PRNG1()*(q-1));
sites2[i][j] |= --spin_value;
num_other_spins++;
}
else if ( PRNG1() < Q_STATE_POTTS_UP_SPIN_PROBABILITY )
// Spin is up.
num_up_spins++;
else
{
// Set spin to another direction.
spin_value = 2 + (int)(PRNG1()*(q-1));
sites2[i][j] |= --spin_value;
num_other_spins++;
}
}
}
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;
// Site is occupied.
sites3[i][j][k] &= ~spin_value_mask; // Set spin to up.
if ( initial_magnetization == 1.0 )
// Spin is up.
num_up_spins++;
else if ( initial_magnetization == -1.0/(q-1))
{
// Set spin to another direction.
spin_value = 2 + (int)(PRNG1()*(q-1));
sites3[i][j][k] |= --spin_value;
num_other_spins++;
}
else if ( PRNG1() < Q_STATE_POTTS_UP_SPIN_PROBABILITY )
// Spin is up.
num_up_spins++;
else
{
// Set spin to another direction.
spin_value = 2 + (int)(PRNG1()*(q-1));
sites3[i][j][k] |= --spin_value;
num_other_spins++;
}
}
}
}
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;
// Site is occupied.
sites4[i][j][k][l] &= ~spin_value_mask; // Set spin to up.
if ( initial_magnetization == 1.0 )
// Spin is up.
num_up_spins++;
else if ( initial_magnetization == -1.0/(q-1))
{
// Set spin to another direction.
spin_value = 2 + (int)(PRNG1()*(q-1));
sites4[i][j][k][l] |= --spin_value;
num_other_spins++;
}
else if ( PRNG1() < Q_STATE_POTTS_UP_SPIN_PROBABILITY )
// Spin is up.
num_up_spins++;
else
{
// Set spin to another direction.
spin_value = 2 + (int)(PRNG1()*(q-1));
sites4[i][j][k][l] |= --spin_value;
num_other_spins++;
}
}
}
}
}
break;
}
if ( adjust_zero_initial_magnetization )
if ( num_other_spins - num_up_spins == num_spins/2 )
{
flip_initial_spin_to_up();
num_up_spins++;
num_other_spins--;
}
abs_magnetization = num_up_spins - num_other_spins;
if ( internal_energy_measured )
// Since the bonds have already been established we can now:
abs_internal_energy = get_abs_internal_energy();
}