Skip to content

Commit a5ebdba

Browse files
dc250601dc250601
committed
update
1 parent 0e529db commit a5ebdba

24 files changed

+3181
-3058
lines changed

‎.ipynb_checkpoints/Untitled-checkpoint.ipynb

+410-599
Large diffs are not rendered by default.

‎Untitled.ipynb

+602-2,320
Large diffs are not rendered by default.

‎hmc_new/control.c

+69
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
// Front code for lattice Gross-Neveu model
2+
// in 2-dim using HMC / Kramer's algorithm
3+
4+
#define CONTROL
5+
#include "include_gn.h"
6+
7+
int main(int argc,char *argv[])
8+
{
9+
10+
int prompt;
11+
int algoflagi,todo;
12+
int i,j,k,l,m,n,nf;
13+
int krflag,numtrial,numaccp;
14+
double avsigma,dtime;
15+
site *s;
16+
17+
FILE *f_sigma;
18+
FILE *fsig,*fpbp;
19+
20+
// Initialize
21+
dtime=-dclock();
22+
prompt=setup();
23+
readin(prompt);
24+
25+
// malloc storages
26+
Nmeas=(int)((KRtraj*KRmax)/Measlen);
27+
av_sigma=(double *)malloc(Nmeas*sizeof(double));
28+
29+
/*****************************************/
30+
// read binary sigma fields, if requested
31+
// else populate it with flat / gaussian RN
32+
if(readflag==FRESH) coldlat();
33+
else if(readflag==FRESHHOT) hotlat();
34+
else{
35+
f_sigma=fopen(read_binary_sigma,"r");
36+
read_bin_double(f_sigma);
37+
fclose(f_sigma);
38+
}
39+
l=numtrial=numaccp=0;
40+
avsigma=0;
41+
42+
for(todo=0;todo<NMdtraj;todo++){
43+
44+
numtrial+=hmc();
45+
numaccp++;
46+
47+
if((todo%Measlen)==0){
48+
printf("trajectory= %d\t av_sigma= %e\n",todo,fabs(average_sigma()));
49+
}
50+
51+
} // todo-loop
52+
53+
printf("\nAcceptance rate = %.4lf\n",(double)numaccp/(double)numtrial);
54+
55+
// save sigma if requested
56+
if(saveflag!=FORGET){
57+
f_sigma=fopen(save_binary_sigma,"w");
58+
save_bin_double(f_sigma);
59+
fclose(f_sigma);
60+
}
61+
62+
// Mark time
63+
dtime+=dclock();
64+
printf("2-dim KS Gross-Neveu completed!\n");
65+
printf("Time = %e seconds\n",dtime);
66+
fflush(stdout);
67+
68+
} // end of main
69+

‎hmc_new/defines.h

+26
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
#ifndef _DEFINES_H
2+
#define _DEFINES_H
3+
4+
#define Xup 0
5+
#define Tup 1
6+
#define Tdn 2
7+
#define Xdn 3
8+
#define Ndirs 4
9+
#define OPPdir(dir) (3-(dir)) // Opposite direction
10+
11+
#define COLDLAT 0
12+
#define HOTLAT 1
13+
#define PLUS 1
14+
#define MINUS -1
15+
16+
#define FRESH 0
17+
#define FORGET 1
18+
#define FRESHHOT 2
19+
#define MAXFILENAME 256
20+
21+
#define HMC 49
22+
#define KMC 50
23+
24+
#define FORALLSITES(i,s) for(i=0,s=lattice;i<Volume;i++,s++)
25+
26+
#endif // _DEFINES_H

‎hmc_new/dslash_gn.c

+147
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
// Dirac matrix multiplication for Gross-Neveu model
2+
// Multiplication involves 4 different scenario :
3+
// M * pseudofermion = pseudofermion => matp2p()
4+
// M * pseudofermion = scalar (double) => matp2d()
5+
// M * scalar (double) = pseudofermion => matd2p()
6+
// M * scalar (double) = scalar (double)=> matd2d()
7+
8+
#include "include_gn.h"
9+
10+
////////////////////////////////////////////////
11+
void matp2p(field_offset src,field_offset dest,int isign,int flav)
12+
{
13+
14+
int i;
15+
site *s;
16+
17+
gather(src,Xup,gen_pt);
18+
FORALLSITES(i,s){
19+
((psferm *)F_PT(s,dest))->f[flav]=(s->sign)*0.5*
20+
((psferm *)gen_pt[i])->f[flav];
21+
}
22+
gather(src,Xdn,gen_pt);
23+
FORALLSITES(i,s){
24+
((psferm *)F_PT(s,dest))->f[flav]-=(s->sign)*0.5*
25+
((psferm *)gen_pt[i])->f[flav];
26+
}
27+
28+
gather(src,Tup,gen_pt);
29+
FORALLSITES(i,s){
30+
((psferm *)F_PT(s,dest))->f[flav]+=0.5*
31+
((psferm *)gen_pt[i])->f[flav];
32+
}
33+
gather(src,Tdn,gen_pt);
34+
FORALLSITES(i,s){
35+
((psferm *)F_PT(s,dest))->f[flav]-=0.5*
36+
((psferm *)gen_pt[i])->f[flav];
37+
}
38+
39+
FORALLSITES(i,s){
40+
((psferm *)F_PT(s,dest))->f[flav]=(isign*
41+
((psferm *)F_PT(s,dest))->f[flav]) +
42+
((s->phi)*((psferm *)F_PT(s,src))->f[flav]);
43+
}
44+
45+
} // matp2p()
46+
47+
////////////////////////////////////////////////
48+
void matp2d(field_offset src,field_offset dest,int isign,int flav)
49+
{
50+
51+
int i;
52+
site *s;
53+
54+
gather(src,Xup,gen_pt);
55+
FORALLSITES(i,s){
56+
*((double *)F_PT(s,dest))=(s->sign)*0.5*
57+
((psferm *)gen_pt[i])->f[flav];
58+
}
59+
gather(src,Xdn,gen_pt);
60+
FORALLSITES(i,s){
61+
*((double *)F_PT(s,dest))-=(s->sign)*0.5*
62+
((psferm *)gen_pt[i])->f[flav];
63+
}
64+
65+
gather(src,Tup,gen_pt);
66+
FORALLSITES(i,s){
67+
*((double *)F_PT(s,dest))+=0.5*((psferm *)gen_pt[i])->f[flav];
68+
}
69+
gather(src,Tdn,gen_pt);
70+
FORALLSITES(i,s){
71+
*((double *)F_PT(s,dest))-=0.5*((psferm *)gen_pt[i])->f[flav];
72+
}
73+
74+
FORALLSITES(i,s){
75+
*((double *)F_PT(s,dest))=(isign*(*((double *)F_PT(s,dest)))) +
76+
((s->phi)*((psferm *)F_PT(s,src))->f[flav]);
77+
}
78+
79+
} // matp2d()
80+
81+
////////////////////////////////////////////////
82+
void matd2p(field_offset src,field_offset dest,int isign,int flav)
83+
{
84+
85+
int i;
86+
site *s;
87+
88+
gather(src,Xup,gen_pt);
89+
FORALLSITES(i,s){
90+
((psferm *)F_PT(s,dest))->f[flav]=(s->sign)*0.5*
91+
(*((double *)gen_pt[i]));
92+
}
93+
gather(src,Xdn,gen_pt);
94+
FORALLSITES(i,s){
95+
((psferm *)F_PT(s,dest))->f[flav]-=(s->sign)*0.5*
96+
(*((double *)gen_pt[i]));
97+
}
98+
99+
gather(src,Tup,gen_pt);
100+
FORALLSITES(i,s){
101+
((psferm *)F_PT(s,dest))->f[flav]+=0.5*(*((double *)gen_pt[i]));
102+
}
103+
gather(src,Tdn,gen_pt);
104+
FORALLSITES(i,s){
105+
((psferm *)F_PT(s,dest))->f[flav]-=0.5*(*((double *)gen_pt[i]));
106+
}
107+
108+
FORALLSITES(i,s){
109+
((psferm *)F_PT(s,dest))->f[flav]=(isign*
110+
((psferm *)F_PT(s,dest))->f[flav]) +
111+
((s->phi)*(*((double *)F_PT(s,src))));
112+
}
113+
114+
} // matd2p()
115+
116+
////////////////////////////////////////////////
117+
void matd2d(field_offset src,field_offset dest,int isign)
118+
{
119+
120+
int i;
121+
site *s;
122+
123+
gather(src,Xup,gen_pt);
124+
FORALLSITES(i,s){
125+
*((double *)F_PT(s,dest))=(s->sign)*0.5*(*((double *)gen_pt[i]));
126+
}
127+
gather(src,Xdn,gen_pt);
128+
FORALLSITES(i,s){
129+
*((double *)F_PT(s,dest))-=(s->sign)*0.5*(*((double *)gen_pt[i]));
130+
}
131+
132+
gather(src,Tup,gen_pt);
133+
FORALLSITES(i,s){
134+
*((double *)F_PT(s,dest))+=0.5*(*((double *)gen_pt[i]));
135+
}
136+
gather(src,Tdn,gen_pt);
137+
FORALLSITES(i,s){
138+
*((double *)F_PT(s,dest))-=0.5*(*((double *)gen_pt[i]));
139+
}
140+
141+
FORALLSITES(i,s){
142+
*((double *)F_PT(s,dest))=(isign*(*((double *)F_PT(s,dest)))) +
143+
((s->phi)*(*((double *)F_PT(s,src))));
144+
}
145+
146+
} // matd2d()
147+

‎hmc_new/hamiltonian.c

+38
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
// Gross-Neveu Hamiltonian for KS fermion
2+
// (eqn. 4.8, mythesis.pdf)
3+
// H = sum_x [ sigma^2/2g^2 + p^2/2 + sum_f zeta_f^d (A^dA)^{-1} zeta_f ]
4+
// hflag == skip inversion, flag = A^{-1} or (A^d A)^{-1}
5+
6+
#include "include_gn.h"
7+
8+
double hamil(int hflag,int flag)
9+
{
10+
11+
int i,nf;
12+
double hml=0,psf_hml;
13+
site *s;
14+
15+
// sigma^2/2g^2 + p^2/2, phi=sigma
16+
FORALLSITES(i,s) hml+=( (0.5/(G*G))*(s->phi)*(s->phi) +
17+
(0.5)*(s->pi)*(s->pi) );
18+
19+
if(hflag!=0){
20+
21+
// sum_f zeta_f^d (A^dA)^{-1} zeta_f
22+
for(nf=0;nf<Nf;nf++){
23+
24+
// eta_f=(A^dagger A)^{-1} zeta_f
25+
congrad(F_OFFSET(zeta),F_OFFSET(eta),CgiterH,ResidueH,flag,nf);
26+
FORALLSITES(i,s) hml+=(s->zeta.f[nf])*(s->eta.f[nf]);
27+
}
28+
} // if-ends
29+
else{
30+
for(nf=0;nf<Nf;nf++){
31+
FORALLSITES(i,s) hml+=(s->zeta.f[nf])*(s->eta.f[nf]);
32+
}
33+
}
34+
35+
return(hml);
36+
37+
} // end of hamil()
38+

‎hmc_new/hmc.c

+70
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
// Kramer equation algorithm --
2+
// 1) starts with stochastic initialization of momentum
3+
// 2) next is momentum initial half step
4+
// 3) full single step for sigma and momentum
5+
// 4) momentum final half step
6+
// 5) Hamiltonian calculation before and after MD step
7+
// followed by Metropolis accept/reject
8+
9+
#include "include_gn.h"
10+
11+
int hmc(void)
12+
{
13+
14+
int i,j,nf,ctr;
15+
double term1,term2,xx,z;
16+
double hold,hnew,deltah;
17+
site *s;
18+
19+
xx=1.0;z=0.0;
20+
for(ctr=0;ctr<(max_hmc_iter && xx>z);ctr++){ // loop till we get
21+
// an acceptance or max_iter
22+
23+
// initialize phi=sigma and momentum=pi
24+
FORALLSITES(i,s){
25+
s->phi=s->sigma;
26+
s->pi=gsl_ran_gaussian(r,1.0);
27+
for(nf=0;nf<Nf;nf++) s->zeta.f[nf]=gsl_ran_gaussian(r,1.0/sqrt(2.0));
28+
}
29+
30+
// initial hamiltonian
31+
hold=hamil(1,1);
32+
//printf("\nOld Hamiltonian = %e\n\n",hold);
33+
34+
// momentum first half-step
35+
piup((double)Step/2.0);
36+
37+
// leapfrog for n-1 steps
38+
for(j=0;j<(Mdstep-1);j++){
39+
FORALLSITES(i,s){ s->phi+=(s->pi*(double)Step);
40+
piup((double)Step);
41+
}
42+
43+
// final n-th step
44+
FORALLSITES(i,s){ s->phi+=(s->pi*(double)Step); }
45+
46+
// momentum final half-step
47+
piup((double)Step/2.0);
48+
49+
// final hamiltonian
50+
hnew=hamil(1,1);
51+
//printf("\nNew Hamiltonian = %e\n\n",hnew);
52+
53+
// Metropolis accept/reject
54+
xx=gsl_ran_flat(r,0.0,1.0);
55+
deltah=hnew-hold;
56+
z=exp(-1.0*deltah);
57+
//printf("\nDeltaH = %e\n",deltah);
58+
59+
// keep looping till xx <= z
60+
61+
} // ctr-loop ends
62+
63+
if(xx<z){ // accepted
64+
//printf("Accepted!\n\n");
65+
FORALLSITES(i,s){ s->sigma=s->phi; }
66+
}
67+
68+
return(ctr);
69+
70+
} // hmc()

0 commit comments

Comments
 (0)