c
===================================================================
program PIGS
c
===================================================================
implicit real*8(a-h,o-z)
parameter (nslicemax=1000)
parameter (nh = 400, xhrange = 6.d00, xhdel = xhrange/nh)
parameter (nsparce = 5)
c
common /var/ x(nslicemax),rho(nslicemax)
common /enesum/ potsum(nslicemax), esum
common /ene/ potenziale(nslicemax), etotale
common /ene2/ potenziale2(nslicemax), etotale2
common /res/ hist(-nh:nh),xnorma,nacc,nall,naccvar
common /spost/ dxmax,dxvar,dxall
common /cost/ cpot,cke
common /cerr/ usupas, usupasmu
common /tempo/ beta, tau
c
dimension irn(4)
c
open (30,file='input.dat',status='old')
open (69,file='Primes',status='old')
read (30,'(4i4)') irn
! seme del generatore di numeri casuali
read (30,*) amass
! massa della particella
read (30,*) beta
! tempo di propagazione
read (30,*) nstep,nbloc,nequil ! numero
di step Monta Carlo
read (30,*) dxmax,dxvar,dxall !
massimo spostamento
read (30,*) nsli,nproiez !
N-tslices tot e necessarie alla proiez
read (69,*) ipr1,ipr2
call setrn(irn(1),irn(2),irn(3),irn(4),ipr1,ipr2,1)
c
c the
program works in atomic units: hbar=1, m_electr=1
c
tau = beta / dble(nsli-1)
alamb = 0.5d00 / amass
! alamb = hbar**2 / 2m
cpot = tau
! costante molt. potential action
cke = 0.25d00 / (alamb*tau) !
costante molt. kinetic action
usupas = 1.d00 / dble(nbloc)
usupasmu = 1.d00 / dble(nbloc-1)
c
xnorma = 0.d00 !
** Zero average value and histogram **
do ih = -nh, nh
hist(ih) = 0.d00
enddo
c
do js=1,nsli
! ** Initial position **
x(js) = 0.d00
! the same position is given to all time slices
potenziale(js) = 0.d00
potenziale2(js) = 0.d00
enddo
etotale = 0.d00
etotale2 = 0.d00
do js=1,nsli-1
rho(js) = densm(cpot,cke,x(js),x(js+1))
enddo
psi1 = funz(x(1))
psi2 = funz(x(nsli))
nacc = 0
naccvar = 0
nall = 0
c
do istep = 1, nequil !
equilibrazione
call metroend(1,psi1,nsli)
call metroend(nsli,psi2,nsli)
call metropigs(nsli,alamb)
call metroall(nsli,psi1,psi2)
enddo
c
nacc = 0
naccvar = 0
nall = 0
do jbloc = 1, nbloc !
RUN
call azero(nsli,iend)
do istep = 1, nstep
if(istep.eq.nstep) iend=1
do isparce = 1, nsparce
call metroend(1,psi1,nsli)
call metroend(nsli,psi2,nsli)
call metropigs(nsli,alamb)
call metroall(nsli,psi1,psi2)
enddo
call estima(iend,nstep,nproiez,nsli,alamb)
enddo
enddo
c
call result(amass,nsli,nstep,nbloc)
call savern(irn(1),irn(2),irn(3),irn(4))
open (96,file='random.out',status='unknown')
write(96,'(4i4)') irn
stop
end
c
c###########
subroutine e funzioni #####################
c
subroutine azero(nsli,iend)
implicit real*8 (a-h,o-z)
parameter (nslicemax=1000)
common /enesum/ potsum(nslicemax), esum
do js=1,nsli
! azzeramento accumulatori
potsum(js) = 0.d00
enddo
esum = 0.d00
iend = 0
return
end
c
subroutine estima(iend,nstep,nproiez,nsli,alamb)
implicit real*8 (a-h,o-z)
parameter (nslicemax=1000)
parameter (nh = 400, xhrange = 6.d00, xhdel = xhrange/nh)
common /var/ x(nslicemax),rho(nslicemax)
common /enesum/ potsum(nslicemax), esum
common /ene/ potenziale(nslicemax), etotale
common /ene2/ potenziale2(nslicemax), etotale2
common /res/ hist(-nh:nh),xnorma,nacc,nall,naccvar
do js=nproiez+1,nsli-nproiez-1
ih = nint(x(js)/xhdel)
if ( abs(ih).le.nh) then
hist(ih) = hist(ih)+1.d00
xnorma = xnorma + 1.d00
endif
enddo
do js=1,nsli
potsum(js) = potsum(js) + field(x(js))
enddo
stimacin= -alamb *(eloc(x(1)) + eloc(x(nsli)))/2.d00
esum=esum+stimacin+(field(x(1))+field(x(nsli)))/2.d00
if(iend.eq.1) then
do js=1,nsli
potmedio =
potsum(js) / dble(nstep)
potenziale(js) = potenziale(js) + potmedio
potenziale2(js)= potenziale2(js)+ potmedio * potmedio
enddo
esum = esum / dble(nstep)
etotale = etotale + esum
etotale2 = etotale2 + esum * esum
endif
return
end
c
subroutine metropigs(nsli,alamb)
implicit real*8 (a-h,o-z)
parameter (nslicemax=1000)
parameter (nh = 400)
common /var/ x(nslicemax),rho(nslicemax)
common /res/ hist(-nh:nh),xnorma,nacc,nall,naccvar
common /spost/ dxmax,dxvar,dxall
common /cost/ cpot,cke
common /tempo/ beta, tau
pi = acos(-1.d00)
dpi=2.d00*pi
do js=2,nsli-1
c set
trial position
xnew = x(js) + (2.d00*rannyu()-1.d00) * dxmax
rsnew = densm(cpot,cke,x(js-1),xnew)
rdnew = densm(cpot,cke,xnew,x(js+1))
c calcolo
della differenza tra l'azione prima e dopo la mossa
delta = rsnew + rdnew - rho(js-1) - rho(js)
c metropolis
per campionare la densita' di probabilita' e^{-S}
if ( dexp(-delta) .gt. rannyu() ) then
c**
... accept **
nacc = nacc + 1
rho(js-1) = densm(cpot,cke,x(js-1),xnew)
rho(js) = densm(cpot,cke,xnew,x(js+1))
x(js) = xnew
endif
enddo
return
end
c
subroutine metroall(nsli,psi1,psi2)
implicit real*8 (a-h,o-z)
parameter (nslicemax=1000)
parameter (nh = 400)
common /var/ x(nslicemax),rho(nslicemax)
common /res/ hist(-nh:nh),xnorma,nacc,nall,naccvar
common /spost/ dxmax,dxvar,dxall
common /cost/ cpot,cke
dimension xold(nslicemax), rhold(nslicemax)
dx = (2.d00*rannyu()-1.d00) * dxall
do js=1,nsli
xold(js) = x(js)
c set
trial position
x(js) = x(js) + dx
enddo
do js=1,nsli-1
rhold(js) = rho(js)
enddo
psi1new = funz(x(1))
psi2new = funz(x(nsli))
delta = 0.d00
do js=1,nsli-1
rho(js) = densm(cpot,cke,x(js),x(js+1))
c calcolo
della differenza tra l'azione prima e dopo la mossa
delta = delta + rho(js) - rhold(js)
c metropolis
per campionare la densita' di probabilita' e^{-S}
enddo
acc = psi1new*psi2new*dexp(-delta)/(psi1*psi2)
if ( acc .gt. rannyu() ) then
c**
... accept **
nall = nall + 1
psi1 = psi1new
psi2 = psi2new
else
c**
... reject and replace **
do js=1,nsli
x(js) = xold(js)
enddo
do js=1,nsli-1
rho(js) = rhold(js)
enddo
endif
c
return
end
c
subroutine metroend(js,psi,nsli)
implicit real*8 (a-h,o-z)
parameter (nslicemax=1000)
parameter (nh = 400)
common /var/ x(nslicemax),rho(nslicemax)
common /res/ hist(-nh:nh),xnorma,nacc,nall,naccvar
common /spost/ dxmax,dxvar,dxall
common /cost/ cpot,cke
if(js.eq.1) then
jsp= 2
rhold = rho(js)
else
jsp=nsli-1
rhold = rho(jsp)
endif
xnew = x(js) + (2.d00*rannyu()-1.d00) * dxvar
psinew = funz(xnew)
rhonew = densm(cpot,cke,xnew,x(jsp))
delta = rhonew - rhold
acc = psinew*dexp(-delta)/psi
c metropolis
if ( acc .gt. rannyu() ) then
c**
... accept **
naccvar = naccvar + 1
x(js) = xnew
psi = psinew
if(js.eq.1) then
rho(js) = rhonew
else
rho(jsp)= rhonew
endif
endif
return
end
c
subroutine result(amass,nsli,nstep,nbloc)
implicit real*8 (a-h,o-z)
parameter (nslicemax=1000)
parameter (nh = 400, xhrange = 6.d00, xhdel = xhrange/nh)
parameter (nsparce = 5)
common /ene/ potenziale(nslicemax), etotale
common /ene2/ potenziale2(nslicemax), etotale2
common /res/ hist(-nh:nh),xnorma,nacc,nall,naccvar
dimension errpot(nslicemax)
do ih = -nh, nh ! ** Normalize average and histogram
**
hist(ih) = hist(ih)/(xnorma*xhdel)
enddo
ratio = dble(nacc) / dble(nsparce*nbloc*nstep*(nsli-2))
ratiovar = dble(naccvar) / dble(nsparce*nbloc*nstep*2)
ratioall = dble(nall) / dble(nsparce*nbloc*nstep)
open(unit=6,file='output.run',status='unknown')
write(6,*) 'Run completed' ! Print run information
write(6,*) 'Acceptance ratio single moves =', ratio
write(6,*) 'Acceptance ratio variational =', ratiovar
write(6,*) 'Acceptance ratio global moves =', ratioall
c**
Print histogram with exact result to compare
**
pi = acos(-1.d00)
alfa = sqrt(amass)
hnorm = sqrt(alfa/pi)
open(unit=3,file='output.psi',status='unknown')
avx = 0.d00
do ih = -nh, nh
xt = dble(ih)*xhdel
write(3,*) xt, hist(ih), hnorm*dexp(-alfa*xt**2)
enddo
close(unit=3)
open(unit=4,file='output.pot',status='unknown')
do i = 1,nsli
poti = potenziale(i)
poti2 = potenziale2(i)
errpot(i) = error(poti,poti2)
potenziale(i) = potenziale(i)/dble(nbloc)
enddo
errtot = error(etotale,etotale2)
etotale = etotale/dble(nbloc)
do i = 1,nsli
write(4,*) i,potenziale(i),errpot(i)
enddo
write(6,*) 'energia totale =',etotale,' +/-',errtot
close(unit=4)
return
end
real*8 function densm(cp,ck,x1,x2)
implicit real*8(a-h,o-z)
common /tempo/ beta, tau
densm = ck*(x1-x2)**2+0.5*cp*(field(x1)+field(x2))
c
fact = ck*tau/dsinh(tau) ! exact harmonic oscillator
density matrix
c
densm = fact*((x1*x1+x2*x2)*dcosh(tau)-2*x1*x2)
return
end
real*8 function funz(x)
implicit real*8(a-h,o-z)
funz = 1.d00/dcosh(x)
c
funz = 1.d00
c
funz = dexp(-0.5d00*x**2)
return
end
real*8 function field(x)
implicit real*8(a-h,o-z)
field = 0.5d00*x**2
return
end
real*8 function eloc(x)
implicit real*8(a-h,o-z)
eloc = -1.d00+2.d00*(dsinh(x)**2)/(dcosh(x)**2)
c
eloc = 0.d00
c
eloc = x*x-1.d00
return
end
real*8 function error(xacm,xacm2)
implicit real*8 (a-h,o-z)
common /cerr/ usupas, usupasmu
xacmsw = xacm * xacm * usupas * usupas
xacm2 = xacm2 * usupas
if (xacm2 .lt. xacmsw ) then
error = 0.d00
else
error = sqrt( (xacm2 - xacmsw) * usupasmu )
endif
return
end
c###################
GENERATORE NUMERI CASUALI ####################
real*8 function rannyu()
implicit real*8 (a-h,o-z)
parameter(ooto12=1.d0/4096.d0)
parameter(itwo12=4096)
common /rnyucm/ l1,l2,l3,l4,n1,n2,n3,n4
data m1,m2,m3,m4 / 502,1521,4071,2107/
i1=l1*m4+l2*m3+l3*m2+l4*m1 + n1
i2=l2*m4+l3*m3+l4*m2 + n2
i3=l3*m4+l4*m3 + n3
i4=l4*m4 + n4
l4=mod(i4,itwo12)
i3=i3+i4/itwo12
l3=mod(i3,itwo12)
i2=i2+i3/itwo12
l2=mod(i2,itwo12)
l1=mod(i1+i2/itwo12,itwo12)
rannyu=ooto12*(dble(l1)+
+ ooto12*(dble(l2)+
+ ooto12*(dble(l3)+
+ ooto12*(dble(l4)))))
return
end
subroutine setrn(j1,j2,j3,j4,jsetac3,jsetac4,jsetac)
implicit real*8 (a-h,o-z)
common /rnyucm/ l1,l2,l3,l4,n1,n2,n3,n4
l1 = j1
l2 = j2
l3 = j3
l4 = j4
n1=0
n2=0
n3=2896
n4=1263
if(jsetac.ne.0)then
n3=jsetac3
n4=jsetac4
endif
return
end
subroutine savern(j1,j2,j3,j4)
implicit real*8 (a-h,o-z)
common /rnyucm/ l1,l2,l3,l4,n1,n2,n3,n4
j1 = l1
j2 = l2
j3 = l3
j4 = l4
return
end