forked from golddoushi/mcsolver
-
Notifications
You must be signed in to change notification settings - Fork 0
/
isingLib.c
183 lines (159 loc) · 6.25 KB
/
isingLib.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
#include "Python.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct Orb
{
int id;
float spin;
int nlink;
float *linkStrength;
int inBlock;
struct Orb **linkedOrb;
}Orb;
void establishLattice(Orb *lattice, int totOrbs, float initSpin[totOrbs], int nlink, float linkStrength[nlink]){
for(int i=0;i<totOrbs;i++){
lattice[i].id=i;
lattice[i].spin=initSpin[i];
lattice[i].nlink=nlink;
lattice[i].linkStrength=linkStrength;
}
}
void establishLinking(Orb *lattice, int totOrbs, int nlink, int linkedOrb[totOrbs][nlink]){
for(int iorb=0;iorb<totOrbs;iorb++){
lattice[iorb].linkedOrb=(Orb**)malloc(nlink*sizeof(Orb*));
for(int ilink=0;ilink<nlink;ilink++){
lattice[iorb].linkedOrb[ilink]=lattice+linkedOrb[iorb][ilink];
}
}
}
float getCorrEnergy(Orb *source){
float corr=0;
for(int i=0;i<source->nlink;i++){
corr+=(source->linkStrength[i])*(source->spin)*(source->linkedOrb[i]->spin);
}
return corr;
}
int expandBlock(int*beginIndex, int*endIndex, Orb *buffer[], int*blockLen, Orb *block[]){
//printf(" Buffer: now start and end pt is %d, %d.\n",*beginIndex, *endIndex);
if(*beginIndex>*endIndex) return 0;
// FIFO
Orb *outOrb=buffer[*beginIndex];
*beginIndex+=1; // pop out the first element
//FILO
//Orb *outOrb=buffer[*endIndex];
//*endIndex-=1; // pop out the last element
int i;
for(i=0;i<outOrb->nlink;i++){
Orb *linkedOrb=outOrb->linkedOrb[i];
//printf(" considering the %d orb which is linking to %d orb, it is %d in block \n", linkedOrb->id, outOrb->id, linkedOrb->inBlock);
if(linkedOrb->inBlock==0){
float corr=(outOrb->linkStrength[i])*(outOrb->spin)*(linkedOrb->spin); // bond strength
//printf(" since it is not in block thus we calc. the correlation energy is %f\n",corr);
if(corr<0 && (1-exp(2*corr))>rand()/32767.0){
//printf(" -->>fortunately it is added to block with possibility %f\n",(1-exp(2*corr)));
// update block
*blockLen+=1;
block[*blockLen-1]=linkedOrb;
linkedOrb->inBlock=1; // register into block
// update buffer
*endIndex+=1;
buffer[*endIndex]=linkedOrb;
}
}
}
return 1;
}
void blockUpdate(int totOrbs, Orb lattice[], float*p_energy, float*p_totSpin){
//printf("one block update step is initializaing...\n");
Orb *block[totOrbs];
Orb *buffer[totOrbs];
int seedID=rand()%totOrbs;
block[0]=lattice+seedID;
buffer[0]=lattice+seedID;
block[0]->inBlock=1;
int beginIndex=0, endIndex=0, blockLen=1, i;
int *p_beginIndex=&beginIndex, *p_endIndex=&endIndex, *p_blockLen=&blockLen;
//printf("the seed Orb is %d\n",block[0]->id);
while (expandBlock(p_beginIndex, p_endIndex, buffer, p_blockLen, block)==1)
{
//printf(" Block size is %d\n",*p_blockLen);
}
for(i=0;i<*p_blockLen;i++){
block[i]->spin*=-1;
block[i]->inBlock=0;
*p_totSpin+=(block[i]->spin*2);
}
*p_energy=0.;
for(i=0;i<totOrbs;i++){
*p_energy+=getCorrEnergy(lattice+i);
}
*p_energy/=2;
}
void localUpdate(int totOrbs, Orb lattice[], float *p_energy, float *p_totSpin){
int seedID=rand()%totOrbs;
float corr=2*getCorrEnergy(lattice+seedID);
if(corr>=0){
lattice[seedID].spin*=-1;
*p_totSpin+=(lattice[seedID].spin*2);
*p_energy-=corr;
}else if (exp(corr)>rand()/32767.0){
lattice[seedID].spin*=-1;
*p_totSpin+=(lattice[seedID].spin*2);
*p_energy-=corr;
}
}
PyObject * blockUpdateMC(int totOrbs, float initSpin[totOrbs], int nthermal, int nsweep,
int nlink, float linkStrength[nlink], int linkedOrb[totOrbs][nlink]){
//printf("hello block algorithm!\n");
// initialize lattice
Orb lattice[totOrbs];
establishLattice(lattice, totOrbs, initSpin, nlink, linkStrength);
establishLinking(lattice, totOrbs, nlink, linkedOrb);
// initialize measurement
float energy=0, totSpin=0;
float *p_energy=&energy, *p_totSpin=&totSpin;
for(int i=0;i<totOrbs;i++) totSpin+=initSpin[i];
// initialize block
for(int i=0;i<totOrbs;i++) lattice[i].inBlock=0;
for(int i=0;i<nthermal;i++) blockUpdate(totOrbs, lattice, p_energy, p_totSpin); //thermalization
// printf("start sweeping\n");
PyObject *spinData, *energyData;
spinData=PyTuple_New(nsweep);
energyData=PyTuple_New(nsweep);
for(int i=0;i<nsweep;i++){
blockUpdate(totOrbs, lattice, p_energy, p_totSpin);
PyTuple_SetItem(spinData, i, PyFloat_FromDouble(*p_totSpin));
PyTuple_SetItem(energyData, i, PyFloat_FromDouble(*p_energy));
}
PyObject *Data;
Data=PyTuple_New(2);
PyTuple_SetItem(Data, 0, spinData);
PyTuple_SetItem(Data, 1, energyData);
return Data;
}
PyObject * localUpdateMC(int totOrbs, float initSpin[totOrbs], int nthermal, int nsweep,
int nlink, float linkStrength[nlink], int linkedOrb[totOrbs][nlink]){
// initialize lattice
Orb lattice[totOrbs];
establishLattice(lattice, totOrbs, initSpin, nlink, linkStrength);
establishLinking(lattice, totOrbs, nlink, linkedOrb);
// initialize measurement
float energy=0, totSpin=0;
float *p_energy=&energy, *p_totSpin=&totSpin;
for(int i=0;i<totOrbs;i++) totSpin+=initSpin[i];
for(int i=0;i<totOrbs*nthermal;i++) localUpdate(totOrbs, lattice, p_energy, p_totSpin); //thermalization
PyObject *spinData, *energyData;
spinData=PyTuple_New(nsweep);
energyData=PyTuple_New(nsweep);
for(int i=0;i<nsweep;i++){
for(int j=0;j<totOrbs;j++) localUpdate(totOrbs, lattice, p_energy, p_totSpin);
PyTuple_SetItem(spinData, i, PyFloat_FromDouble(*p_totSpin));
PyTuple_SetItem(energyData, i, PyFloat_FromDouble(*p_energy));
}
PyObject *Data;
Data=PyTuple_New(2);
PyTuple_SetItem(Data, 0, spinData);
PyTuple_SetItem(Data, 1, energyData);
return Data;
}