c
===================================================================
program PIGS
c
===================================================================
implicit real*8(a-h,o-z)
parameter (nslicemax=1000)
c
common /var/
x(3,nslicemax),rho(nslicemax)
common /enesum/ potsum(nslicemax), esum
common /ene/
potenziale(nslicemax), etotale
common /ene2/
potenziale2(nslicemax), etotale2
common /res/
nacc,nall,naccvar
common /spost/ dxmax,dxvar,dxall
common /cost/ cpot,cke
common /cerr/ usupas,
usupasmu
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
do
js=1,nsli
! ** Initial position **
x(1,js) = 0.d00 ! the
same position is given to all time slices
x(2,js) = 0.d00
x(3,js) = 0.d00
potenziale(js) = 0.d00
potenziale2(js) = 0.d00
enddo
etotale = 0.d00
etotale2 = 0.d00
do js=1,nsli-1
x1 = x(1,js)
y1 = x(2,js)
z1 = x(3,js)
x2 = x(1,js+1)
y2 = x(2,js+1)
z2 = x(3,js+1)
rho(js) =
densm(cpot,cke,x1,y1,z1,x2,y2,z2)
enddo
x1 = x(1,1)
y1 = x(2,1)
z1 = x(3,1)
psi1 = funz(x1,y1,z1)
x2 = x(1,nsli)
y2 = x(2,nsli)
z2 = x(3,nsli)
psi2 = funz(x2,y2,z2)
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
call
metroend(1,psi1,nsli)
call
metroend(nsli,psi2,nsli)
call
metropigs(nsli,alamb)
call
metroall(nsli,psi1,psi2)
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)
common /var/
x(3,nslicemax),rho(nslicemax)
common /enesum/ potsum(nslicemax), esum
common /ene/
potenziale(nslicemax), etotale
common /ene2/
potenziale2(nslicemax), etotale2
common /res/
nacc,nall,naccvar
do js=1,nsli
x1 = x(1,js)
y1 = x(2,js)
z1 = x(3,js)
potsum(js) = potsum(js) +
field(x1,y1,z1)
enddo
x1 = x(1,1)
y1 = x(2,1)
z1 = x(3,1)
x2 = x(1,nsli)
y2 = x(2,nsli)
z2 = x(3,nsli)
stimacin= -alamb
*(eloc(x1,y1,z1)+eloc(x2,y2,z2))/2.d00
esum=esum+stimacin+(field(x1,y1,z1)+field(x2,y2,z2))/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)
common /var/
x(3,nslicemax),rho(nslicemax)
common /res/
nacc,nall,naccvar
common /spost/ dxmax,dxvar,dxall
common /cost/ cpot,cke
pi = acos(-1.d00)
dpi=2.d00*pi
do js=2,nsli-1
c set trial position
xnew = x(1,js) +
(2.d00*rannyu()-1.d00) * dxmax
ynew = x(2,js) +
(2.d00*rannyu()-1.d00) * dxmax
znew = x(3,js) +
(2.d00*rannyu()-1.d00) * dxmax
x1 =
x(1,js-1)
y1 =
x(2,js-1)
z1 =
x(3,js-1)
rsnew =
densm(cpot,cke,x1,y1,z1,xnew,ynew,znew)
x1 =
x(1,js+1)
y1 =
x(2,js+1)
z1 =
x(3,js+1)
rdnew =
densm(cpot,cke,xnew,ynew,znew,x1,y1,z1)
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) = rsnew
rho(js) = rdnew
x(1,js) = xnew
x(2,js) = ynew
x(3,js) = znew
endif
enddo
return
end
c
subroutine metroall(nsli,psi1,psi2)
implicit real*8 (a-h,o-z)
parameter (nslicemax=1000)
common /var/
x(3,nslicemax),rho(nslicemax)
common /res/
nacc,nall,naccvar
common /spost/ dxmax,dxvar,dxall
common /cost/ cpot,cke
dimension xold(3,nslicemax),
rhold(nslicemax)
dx = (2.d00*rannyu()-1.d00) * dxall
dy = (2.d00*rannyu()-1.d00) * dxall
dz = (2.d00*rannyu()-1.d00) * dxall
do js=1,nsli
xold(1,js) = x(1,js)
xold(2,js) = x(2,js)
xold(3,js) = x(3,js)
c set trial position
x(1,js) = x(1,js) + dx
x(2,js) = x(2,js) + dy
x(3,js) = x(3,js) + dz
enddo
do js=1,nsli-1
rhold(js) = rho(js)
enddo
psi1new = funz(x(1,1),x(2,1),x(3,1))
psi2new =
funz(x(1,nsli),x(2,nsli),x(3,nsli))
delta = 0.d00
do js=1,nsli-1
x1 = x(1,js)
y1 = x(2,js)
z1 = x(3,js)
x2 = x(1,js+1)
y2 = x(2,js+1)
z2 = x(3,js+1)
rho(js) =
densm(cpot,cke,x1,y1,z1,x2,y2,z2)
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(1,js) =
xold(1,js)
x(2,js) =
xold(2,js)
x(3,js) =
xold(3,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)
common /var/
x(3,nslicemax),rho(nslicemax)
common /res/
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(1,js) +
(2.d00*rannyu()-1.d00) * dxvar
ynew = x(2,js) +
(2.d00*rannyu()-1.d00) * dxvar
znew = x(3,js) +
(2.d00*rannyu()-1.d00) * dxvar
psinew = funz(xnew,ynew,znew)
x1 = x(1,jsp)
y1 = x(2,jsp)
z1 = x(3,jsp)
rhonew =
densm(cpot,cke,xnew,ynew,znew,x1,y1,z1)
delta = rhonew - rhold
acc = psinew*dexp(-delta)/psi
c metropolis
if ( acc .gt. rannyu() ) then
c** ... accept **
naccvar = naccvar + 1
x(1,js) = xnew
x(2,js) = ynew
x(3,js) = znew
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)
common /ene/
potenziale(nslicemax), etotale
common /ene2/ potenziale2(nslicemax),
etotale2
common /res/ nacc,nall,naccvar
dimension errpot(nslicemax)
ratio = dble(nacc) /
dble(nbloc*nstep*(nsli-2))
ratiovar = dble(naccvar) /
dble(nbloc*nstep*2)
ratioall = dble(nall) / dble(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
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,y1,z1,x2,y2,z2)
implicit real*8(a-h,o-z)
dmkin =
ck*((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
dmpot =
0.5d0*cp*(field(x1,y1,z1)+field(x2,y2,z2))
densm = dmkin + dmpot
return
end
real*8 function funz(x,y,z)
implicit real*8(a-h,o-z)
funz = 1.d0
return
end
real*8 function field(x,y,z)
implicit real*8(a-h,o-z)
field = 0.5d0*(x**2 + y**2 + z**2)
return
end
real*8 function eloc(x,y,z)
implicit real*8(a-h,o-z)
eloc = 0.d0
return
end
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