#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#define SIZE 512
/*Num of sites*/
#define Nx 24
/*Num of sites in x direction*/
#define Ny 24
/*Num of sites in y direction*/
#define Nz 24
/*Num of sites in z direction*/
#define num_clust 5
/*number of total cluster*/
double hamilccu[SIZE][SIZE]; /*charge carrier hamil calculated in matrixcc (spin up)*/
double hamilccd[SIZE][SIZE]; /*charge carrier hamil calculated in matrixcc (spin
down)*/
double jij[SIZE][SIZE]; /*matrix w/overlaps calculated in jfind*/
double hopp[SIZE][SIZE]; /*the coefs of the hopping matrix (jfind)*/
double sofib[SIZE]; /*spin at site i*/
double sofia[SIZE]; /*spin at site i*/
double fdisu[SIZE]; /*fermi dist. for spin ups (func. of ener)*
/
double fdisd[SIZE]; /*fermi dist. for spin downs (func. of ener
double enerup[SIZE];
/*spin up eigenvalues calculated from eigen*/
double enerdown[SIZE]; /*spin down eigenvalues calculated from eigen*/
int
numpart[num_clust+1];
/*number of particles in cluster*/
void
interact(double x[SIZE], double y[SIZE], double z[SIZE]); /*findsdistance between sites*/
void eigen();
/*findseigenvalues and vectors of Hamil*/
void
dsyevd_(); /*NAG
routine used in eigen*/
void
brillouin(double x, double S, int site1); /*<S(mn)>(T)*/
void jfind(double
radial[SIZE][SIZE]); /*findsJij and hopp*/
void
mnplace(); /*puts Mn
atoms on lattice*/
double weight(int x);
double chempf(double ywant, double beta); /*Finds chem pot
based on #holes*/
double nhole(double chempot, double beta); /*Uses fer
mi dist. to find # of holes*/
double nint(double x); /*rounds x
to nearest integer*/
int mncheck(double
error); /*check
s error bounds*/
int
prob_dens(void);
int rand2(double
a, double b, int c);
int test(int n[SIZE],int nc);
int twopart(int
atm_num, double x[SIZE], double y[SIZE], double z[SIZE]);
int fcc(int a, int
b, int c);
FILE *stream;
int norm = 32767;
void saveclus(double x[SIZE], double y[SIZE], double
z[SIZE])
{
char filename[20];
int i, i1, i2,
icluster;
i1 = 0;
i2 = numpart[1];
for(icluster=0;
icluster<num_clust; icluster++)
{
sprintf(filename, "clustero%d.txt", icluster);
stream =
fopen(filename, "w+");
for(i=i1;
i<i2; i++)
fprintf(stream, "%f %f %f\n", x[i], y[i], z[i]);
fclose(stream);
i1 = i2;
i2+=
numpart[icluster+1];
}
return;
}
void savedata(double thm2sp, double asmn, double ascc,
double asccc[num_clust+1]
, double ascmn[num_clust+1], double achd[num_clust+1], int
count2)
{
char filename1[20];
char filename2[20];
char filename3[20];
int i, icluster;
if(count2==0){
stream =
fopen("mnsdata.txt", "w+");
fprintf(stream,
"%f %f\n", thm2sp, asmn);
fclose(stream);
stream =
fopen("ccsdata.txt", "w+");
fprintf(stream,
"%f %f\n", thm2sp, ascc);
fclose(stream);
for(icluster=1;
icluster<=num_clust; icluster++)
{
sprintf(filename1, "asccc%d.txt", icluster);
stream =
fopen(filename1, "w+");
fprintf(stream, "%f %f\n", thm2sp, asccc[icluster]);
fclose(stream);
}
for(icluster=1;
icluster<=num_clust; icluster++)
{
sprintf(filename2, "ascmn%d.txt", icluster);
stream =
fopen(filename2, "w+");
fprintf(stream, "%f %f\n", thm2sp, ascmn[icluster]);
fclose(stream);
}
for(icluster=1;
icluster<=num_clust; icluster++)
{
sprintf(filename3, "achd%d.txt", icluster);
stream =
fopen(filename3, "w+");
fprintf(stream, "%f %f\n", thm2sp, achd[icluster]);
fclose(stream);
}
}
else{
stream =
fopen("mnsdata.txt", "a+");
fprintf(stream,
"%f %f\n", thm2sp, asmn);
fclose(stream);
stream =
fopen("ccsdata.txt", "a+");
fprintf(stream,
"%f %f\n", thm2sp, ascc);
fclose(stream);
for(icluster=1;
icluster<=num_clust; icluster++)
{
sprintf(filename1, "asccc%d.txt", icluster);
stream =
fopen(filename1, "a+");
fprintf(stream, "%f %f\n", thm2sp, asccc[icluster]);
fclose(stream);
}
for(icluster=1;
icluster<=num_clust; icluster++)
{
sprintf(filename2, "ascmn%d.txt", icluster);
stream =
fopen(filename2, "a+");
fprintf(stream, "%f %f\n", thm2sp, ascmn[icluster]);
fclose(stream);
}
for(icluster=1;
icluster<=num_clust; icluster++)
{
sprintf(filename3, "achd%d.txt", icluster);
stream =
fopen(filename3, "a+");
fprintf(stream,
"%f %f\n", thm2sp, achd[icluster]);
fclose(stream);
}
}
return;
}
void mnplace()
{
double x[SIZE],
y[SIZE], z[SIZE];
int atm_num = 0, i,
j, prime;
int atb;
numpart[0]=0;
prime = 8;
x[prime] = 0;
srand(time(0));
do
for(i=1;
i<=num_clust; i++)
numpart[i] =
rand2(1, SIZE/num_clust*2+1, 0);
while(test(numpart,
num_clust));
for(i=1;
i<=num_clust; i++)
for(j=0;
j<numpart[i]; j++)
if(j==0){
do{
x[atm_num]=rand2(0,2*Nx-1,0);
y[atm_num]=rand2(0,2*Ny-1,0);
z[atm_num]=rand2(0,2*Nz-1,0);
atb =
fcc(x[atm_num],y[atm_num],z[atm_num]);
}
prime=atm_num;
atm_num++;
}
else{
do{
x[atm_num]=x[prime]+rand2(-1,1,0)*prob_dens();
y[atm_num]=y[prime]+rand2(-1,1,0)*prob_dens();
z[atm_num]=z[prime]+rand2(-1,1,0)*prob_dens();
atb =
fcc(x[atm_num],y[atm_num],z[atm_num]);
}
while((twopart(atm_num,x,y,z))||(atb));
atm_num++;
}
prime = 0;
for(atm_num=0;
atm_num<SIZE; atm_num++){
x[atm_num]=x[atm_num]/2.;
y[atm_num]=y[atm_num]/2.;
z[atm_num]=z[atm_num]/2.;
}
if(num_clust<15)
saveclus(x, y,
z);
interact(x, y, z);
return;
}
/*makes sure that there are not two particles on the same
site*/
int twopart(int atm_num, double x[SIZE], double y[SIZE],
double z[SIZE])
{
int i,j;
for(i=0;
i<atm_num; i++)
if((x[i]==x[atm_num])&&(y[i]==y[atm_num])&&(z[i]==z[atm_num])){
return 1;
}
return 0;
}
/*ensures particle is on an avaliable fcc site*/
int fcc(int a, int b, int c)
{
int d, e, f, g;
d =
((a%2==0)&&(b%2==0)&&(c%2==0));
e =
((a%2==1)&&(b%2==1)&&(c%2==0));
f =
((a%2==1)&&(b%2==0)&&(c%2==1));
g =
((a%2==0)&&(b%2==1)&&(c%2==1));
if((a<0)||(b<0)||(c<0)||(a>2*Nx)||(b>2*Ny)||(c>2*Nz))
return 1;
if((!d)&&(!e)&&(!f)&&(!g))
return 1;
return 0;
}
/*Makes sure that the total number of particles in the
system is size*/
int test(int n[SIZE], int nc)
{
int i, j;
int tot_part = 0;
for(i=1; i<=nc;
i++)
tot_part +=
n[i];
if(tot_part ==
SIZE)
return 0;
return 1;
}
int prob_dens(void)
{
int y, i;
y = rand();
for(i=0; i<Nx;
i++)
if
(y<norm*weight(i))
return i+1;
return Nx;
}
double weight(int x)
{
double y;
if(x==0)
return .5;
y =
cosh(x)/sinh(x)- .5*cosh(.5*x)/sinh(.5*x)+.5;
return y;
}
int rand2(double a, double b, int c)
{
double x;
double y;
int top, fl;
if(c==1)
srand(time(0));
x = rand();
y = x/norm*(b-a)+a;
y = nint(y);
return y;
}
main()
{
/*dummies*/
int site1; /*dummy variable over
sites*/
int site2; /*dummy variable over
sites*/
int count=1;
int count2=0; /*counts the number of
iterations*/
int flag;
/*constants*/
double gmn =
2.0; /*gyromag. ratio of
Mn*/
double bm =
5.789e-5; /*bohr mag.*/
double j =
1.5e-2; /*overlap energy*/
double spinmn
= 2.5; /*quantum numbers Sz of Mn*/
double boltz =
8.614e-5; /*boltzman constant */
/*parameters*/
double hfld =
0.0; /*magnetic field*/
double tmp = 10.0; /*temperature*/
double hl2mn =
.1; /*ratio of holes to
number of Mn spins (p in paper
)*/
int site_num;
double nh1 = hl2mn*SIZE; /*number of holes (desired)*/
double tnot = 0.01*j/boltz; /*temperature for initial iteration*/
/*unknowns*/
double
aspinmnb; /*average Mn spin
before*/
double
aspinmna; /*average Mn spin
after-- before and after cases should be the same*/
double
aspincc; /*average cc spin*/
double
chempot; /*chemical
potential*/
double nh2; /*number of holes after nhole
runs (calculated)*/
double onsite; /*Onsite energy (epsilon in
paper)*/
double
intener; /*Total Internal
energy*/
/*non fundemental
parameters*/
double beta; /*beta*/
double thm2sp; /*initial thermal energy to spin
energy ratio*/
/*arrays*/
double
effh[SIZE]; /*effective magn.
field at site i*/
double
spincc[SIZE]; /*spin at site i*/
double
clustav[num_clust+1]; /*av. Mn spin
over a cluster*/
double
clustavcc[num_clust+1]; /*av. CC spin
over a cluster*/
double
holedens[SIZE];
double
clusthd[num_clust+1];
/*allocates memory
for the largest array (work) used in eigen*/
double *Work;
Work = (double
*)malloc((8*SIZE*SIZE)*sizeof(double));
aspinmna =
2.5; /*initial aspinmna
guess*/
/*runs mnplace
which finds the distances between sites based on concentration=
>interact*/
mnplace();
for(site1=0;
site1<SIZE; site1++)
sofia[site1]=aspinmna;
for(thm2sp =
tnot*boltz/j ;thm2sp<2; thm2sp+=.02)
{
tmp =
thm2sp*j/boltz;
do
{
if((count>10)&&(aspinmna<1)&&(fabs(aspinmnb-aspinmna)>.05))
{
for(site1=0;
site1<SIZE; site1++)
sofia[site1]=sofia[site1]-(sofib[site1]-sofia[site1]);
}
else
{
aspinmnb =
aspinmna;
for(site1=0; site1<SIZE; site1++)
sofib[site1]=sofia[site1];
}
/*since the
original matrix was lost after eigen ran*/
for(site1=0;
site1<SIZE; site1++)
{
onsite = -gmn*bm*hfld/2.;
for(site2=0; site2<SIZE; site2++)
{
hamilccd[site1][site2]=hopp[site1][site2];
hamilccu[site1][site2]=hopp[site1][site2];
onsite
+= jij[site1][site2]*sofib[site2]/2.;
}
hamilccu[site1][site1]=onsite;
hamilccd[site1][site1]=-onsite;
}
eigen(Work);
beta = 1/(tmp*boltz);
chempot = chempf(nh1,beta);
nh2 = nhole(chempot,beta);
aspincc = 0;
aspinmna = 0;
intener = 0;
/*equation 13
in paper*/
for(site1=0;
site1<SIZE; site1++)
{
spincc[site1]=0;
holedens[site1]=0;
for(site2=0; site2<SIZE; site2++){
spincc[site1]+=(pow(hamilccu[site2][site1],2)*fdisu[site2]-pow(hamil
ccd[site2][site1],2)*fdisd[site2])*.5;
holedens[site1]+=pow(hamilccu[site2][site1],2)*fdisu[site2]+pow(hami
lccd[site2][site1],2)*fdisd[site2];
}
if(holedens[site1]>1)
printf("%d => %f", site1, holedens[site1]);
intener+=enerup[site1]*fdisu[site1]+enerdown[site1]*fdisd[site1];
aspincc+=spincc[site1]/nh2;
}
/*eqaution 5 in
paper defines effective mag field*/
for(site1=0;
site1<SIZE; site1++)
{
effh[site1]=gmn*bm*hfld;
for(site2=0; site2<SIZE; site2++)
effh[site1]-=jij[site1][site2]*spincc[site2];
brillouin(effh[site1]*beta, spinmn, site1);
aspinmna+=sofia[site1]/SIZE;
}
/*
printf("%3d) sb=%f sa=%f h1=%f h2= %f scc=%f T=%f \n", count,
aspinmnb,
aspinmna, nh1, nh2,
aspincc, tmp);*/
count++;
}
while
(mncheck(aspinmna));
/*a check to make
sure that the number of holes I'm getting is in the ball par
k*/
if(fabs(nh1-nh2)>0.01)
printf("nh1
= %f nh2= %f which is unacceptable", nh1, nh2);
if(tmp ==
tnot){
printf("
%d %d ", site1, numpart[site1]);
printf("\n");
}
/*resets count for
any for loops*/
printf("%4.2f
%4.2f", aspinmna, intener);
count = 1;
site_num = 0;
for(site1=1; site1<=num_clust;
site1++){
clustav[site1]=0;
clustavcc[site1]=0;
clusthd[site1]=0;
for(site2=0;
site2<numpart[site1]; site2++)
{
clustavcc[site1]+=spincc[site_num];
clustav[site1]+=sofia[site_num]/numpart[site1];
clusthd[site1]+= holedens[site_num];
site_num++;
}
clustavcc[site1]=clustavcc[site1]/clusthd[site1];
printf("
%4.2f/%4.2f/%4.2f ", clustav[site1], clustavcc[site1], clusthd[site
1]);
}
printf("\n");
savedata(thm2sp,
aspinmna, aspincc, clustavcc, clustav, clusthd, count2);
count2++;
}
return 0;
}
/*def. of brill funct. from stat. mech*/
void brillouin(double x, double S, int site1)
{
/*Defines brill. at
limiting cases*/
if(x < 1e-7)
sofia[site1]=0;
if(x>197)
sofia[site1]=S;
else
sofia[site1] =
(S+.5)*cosh((S+.5)*x)/sinh((S+.5)*x)-.5*cosh(.5*x)/sinh(.5*x)
;
return;
}
/*calculates distances between particle i and j*/
void interact(double x[SIZE], double y[SIZE], double
z[SIZE])
{
double radial[SIZE][SIZE]; /*computed distances between site1 and
site2
from interact.*/
double sx, sy, sz,
sxy, sxz, syz, sxyz, s0;
double xdiff,
ydiff, zdiff;
int s1, s2;
for(s1=0;
s1<SIZE; s1++){
for(s2=0;
s2<SIZE; s2++)
{
ydiff =
fabs(y[s1]-y[s2]);
zdiff =
fabs(z[s1]-z[s2]);
if(xdiff>0.5*Nx)
xdiff = Nx
- xdiff;
if(ydiff>0.5*Ny)
ydiff = Ny
- ydiff;
if(zdiff>0.5*Nz)
zdiff = Nz
- zdiff;
radial[s1][s2]=sqrt(xdiff*xdiff+ydiff*ydiff+zdiff*zdiff);
}
}
jfind(radial);
printf("\n");
return;
}
/*finds overlap matrix and its contractions as well as
hopping matrix*/
void jfind(double radial[SIZE][SIZE])
{
double j =
15e-3;
/*overlap max*/
int site1, site2;
double a =
5.65e-10;
/*a for GaAs*/
double abohr =
7.8e-10;
/*bohr radius*/
double aratio =
a/abohr;
double Ry =
112.4e-3; /*One Rydberg (binding ene
rgy)*/
for(site1=0;
site1<SIZE; site1++)
for(site2=0;
site2<SIZE; site2++)
{
jij[site1][site2] = j*exp(-2*aratio*radial[site1][site2]);
hopp[site1][site2] = 2*(1+aratio*(radial[site1][site2]))*exp(-aratio*rad
ial[site1][site2])*Ry;
}
return;
}
/*calculates eigenvalues and eigenvectors*/
void eigen(double *Work)
{
int site1;
int site2;
int LWORK =
8*SIZE*SIZE;
int N2 = SIZE;
int LDA = SIZE;
int INFO;
char JOB = 'V';
char UPLO = 'L';
int
IWORK[5*SIZE+3];
int LIWORK =
5*N2+3;
/*calculates
eigenvalues for spin up and spin down hamil. respectively*/
dsyevd_(&JOB,
&UPLO, &N2, &hamilccu[0][0], &LDA, &enerup[0],
&Work[0], &LWORK,
&IWORK[0],
&LIWORK, &INFO);
dsyevd_(&JOB,
&UPLO, &N2, &hamilccd[0][0], &LDA, &enerdown[0],
&Work[0], &LWOR
K, &IWORK[0], &LIWORK, &INFO);
return;
}
/*function finds the chem. potential to give the desired
hole to spin ratio*/
double chempf(double ywant, double beta)
{
int count2 = 0;
double lowb =
-10; /*initial lower bound of
argument*/
double uppb =
10; /*initial upper bound of
argument*/
double y;
double x;
do
{
x = (lowb +
uppb)/2.; /*sets arg.= aver. of uppb
and lowb*/
y = nhole(x,
beta); /*defines value of y*/
if(y<ywant) /*if
my value is too large*/
lowb = x;
else
if(y>ywant) /*if my value is
too small*/
uppb =
x;
count2++;
}
while(fabs(lowb - uppb)>0.00000001);
return x;
}
/*figures out how many holes are present for a given
chempot. and temp*/
double nhole(double chempot, double beta)
{
double nh;
int site1, site2;
nh=0;
for(site1=0;
site1<SIZE; site1++)
{
fdisu[site1]=1/(exp(beta*(enerup[site1]-chempot))+1);
fdisd[site1]=1/(exp(beta*(enerdown[site1]-chempot))+1);
nh+=fdisu[site1]+fdisd[site1];
}
return nh;
}
/*figures out whether another iteration is necessary*/
int mncheck(double error)
{
int site1, site2;
int test;
if (error<1)
error = 6e-4;
else
error = 6e-5;
for(site1=0;
site1<SIZE; site1++)
if
(fabs(sofia[site1]-sofib[site1])>error)
return 1;
return 0;
}
/*rounds integer x*/
double nint(double x)
{
double y;
double middleg;
middleg =
(ceil(x)+floor(x))/2.;
if(middleg>x)
y = floor(x);
if(middleg<=x)
y = ceil(x);
return y;
}