mentecaro.for

用fortran写的蒙特卡罗法写的程序 我从国外一个博士论文中下载的 可惜还有一个小错误 希望大家能解决

	else
	if(zi>(W-dpi/2.0)) then
	ovr=1.0
	goto 400
	else
	do k=1,nw
	rxi=xi-xw(k)
	ryi=yi-yw(k)
	rzi=zi-zw(k)
	rxi=rxi-L*anint(rxi/L)
	ryi=ryi-L*anint(ryi/L)
	dri2=rxi*rxi+ryi*ryi+rzi*rzi
	if(dri2>=rminpwsq) then
	dri=sqrt(dri2)
	vpw=vpw+(qpi*qw*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 400
	end if
	end do
	end if
	end if
	! Interaction of positive ions with wall charged sheets
	zdis=abs(zi)
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vpws=nw*(qpi*qw*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	end do ! Close loop over all positive ions
	write(*,*)'No overlap of positive ions'
	write(*,fmt=410)ovr
	! Sum over all negative ions
	do i=1,nn
	xi=xn(i)
	yi=yn(i)
	zi=zn(i)
	! Interaction of negative ions with negative ions
	if(i<nn) then
	do k=i+1,nn
	rxi=xi-xn(k)
	ryi=yi-yn(k)
	rzi=zi-zn(k)
	rxi=rxi-L*anint(rxi/L)
	ryi=ryi-L*anint(ryi/L)
	dri2=rxi*rxi+ryi*ryi+rzi*rzi

	if(dri2>=rminnnsq) then
	dri=sqrt(dri2)
	vnn=vnn+(qn*qn*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 400
	end if
	end do
	end if
	! Interaction of negative ions with positively charged sheets
	do k=1,np
	zdis=abs(zi-zp(k))
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vnps=vnps+(qp(k)*qn*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	end do
	! Interaction of negative ions with negatively charged sheets
	do k=1,nn
	zdis=abs(zi-zn(k))
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vnns=vnns+(qn*qn*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	end do
	! Interaction of negative ions with the walls
	if(zi<(dn/2.0)) then
	ovr=1.0
	goto 400
	else
	if(zi>(W-dn/2.0)) then
	ovr=1.0
	goto 400
	else
	do k=1,nw
	rxi=xi-xw(k)
	ryi=yi-yw(k)
	rzi=zi-zw(k)
	rxi=rxi-L*anint(rxi/L)
	ryi=ryi-L*anint(ryi/L)
	dri2=rxi*rxi+ryi*ryi+rzi*rzi
	if(dri2>=rminnwsq) then
	dri=sqrt(dri2)
	vnw=vnw+(qn*qw*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 400
	end if
	end do
	end if
	end if

	! Interaction of negative ions with wall charged sheets
	zdis=abs(zi)
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vnws=nw*(qn*qw*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	end do ! Close loop over all negative ions
	write(*,*)'No overlap of negative ions'
	write(*,fmt=410)ovr
	vint=vpp+vpps+vpn+vpns+vpw+vpws+vnn+vnps+vnns+vnw+vnws
	write(*,fmt=420)vint
	410 format(1x,f10.6)
	420 format(1x,E14.6)
	400 RETURN
	END SUBROUTINE sumup
	! -----------------------------------------------------------------------------
	SUBROUTINE energyp(dxi,dyi,dzi,ind,vinti,ovr)
	USE mixdat
	IMPLICIT NONE
	integer,intent(in)::ind
	real,intent(in)::dxi,dyi,dzi
	real(kind=8),intent(out)::vinti
	real,intent(out)::ovr
	integer::j,k
	real(kind=8)::vp,vps,vn,vns,vw,vws
	real::xi,yi,zi,zdis,vm
	real::r1,r1aux,r2,r2aux,pos,arctan
	real::EW,aux,fi
	real::rxi,ryi,rzi,dri2,dri
	real::qpi,dpi
	real::rminppsq,rminpnsq,rminnnsq,rminpwsq,rminnwsq
	! Initialization of accumulators
	vinti=0.0
	vp=0.0
	vps=0.0
	vn=0.0
	vns=0.0
	vw=0.0
	vws=0.0
	ovr=0.0
	vm=0.0
	xi=dxi
	yi=dyi
	zi=dzi
	qpi=qp(ind)
	dpi=dp(ind)
	rminpnsq=(dpi/2.0+dn/2.0)*(dpi/2.0+dn/2.0)
	rminpwsq=(dpi/2.0+dw/2.0)*(dpi/2.0+dw/2.0)

	! Interaction with positive ions
	do k=1,np
	if(k/=ind) then
	rxi=xi-xp(k)
	ryi=yi-yp(k)
	rzi=zi-zp(k)
	rxi=rxi-L*anint(rxi/L)
	ryi=ryi-L*anint(ryi/L)
	rminppsq=(dpi/2.0+dp(k)/2.0)*(dpi/2.0+dp(k)/2.0)
	dri2=rxi*rxi+ryi*ryi+rzi*rzi
	if(dri2>=rminppsq) then
	dri=sqrt(dri2)
	vp=vp+(qpi*qp(k)*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 500
	end if
	end if
	end do
	! Interaction with positively charged sheets
	do k=1,np
	zdis=abs(zi-zp(k))
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vps=vps+(qpi*qp(k)*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	end do
	! Interaction with negative ions
	do k=1,nn
	rxi=xi-xn(k)
	ryi=yi-yn(k)
	rzi=zi-zn(k)
	rxi=rxi-L*anint(rxi/L)
	ryi=ryi-L*anint(ryi/L)
	dri2=rxi*rxi+ryi*ryi+rzi*rzi
	if(dri2>=rminpnsq) then
	dri=sqrt(dri2)
	vn=vn+(qpi*qn*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 500
	end if
	end do
	! Interaction with negatively charged sheets
	do k=1,nn
	zdis=abs(zi-zn(k))
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vns=vns+(qpi*qn*elec*elec*fi)/(4.0*pi*eps*epso*L*L)

	end do
	! Interaction with the walls
	if(zi<(dpi/2.0)) then
	ovr=1.0
	goto 500
	else
	if(zi>(W-dpi/2.0)) then
	ovr=1.0
	goto 500
	else
	do k=1,nw
	rxi=xi-xw(k)
	ryi=yi-yw(k)
	rzi=zi-zw(k)
	rxi=rxi-L*anint(rxi/L)
	ryi=ryi-L*anint(ryi/L)
	dri2=rxi*rxi+ryi*ryi+rzi*rzi
	if(dri2>=rminpwsq) then
	dri=sqrt(dri2)
	vw=vw+(qpi*qw*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 500
	end if
	end do
	end if
	end if
	! Interaction with wall charged sheets
	zdis=abs(zi)
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vws=nw*(qpi*qw*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	vinti=vp+vps+vn+vns+vw+vws
	500 RETURN
	END SUBROUTINE energyp
	! --------------------------------------------------------------------------
	SUBROUTINE energyneg(lxi,lyi,lzi,ind,vinti,ovr)
	USE mixdat
	IMPLICIT NONE
	integer,intent(in)::ind
	real,intent(in)::lxi,lyi,lzi
	real(kind=8),intent(out)::vinti
	real,intent(out)::ovr
	integer::j,k
	real::xi,yi,zi
	real(kind=8)::vp,vps,vn,vns,vw,vws
	real::r1aux,r1,r2aux,r2,aux,fi,EW,arctan,pos
	real::rxi,ryi,rzi,zdis,dri,dri2
	real::rminppsq,rminnnsq,rminpnsq,rminpwsq,rminnwsq
	real vm
	! Initialization of accumulators
	vinti=0.0
	vp=0.0
	vps=0.0
	vn=0.0
	vns=0.0
	vw=0.0
	vws=0.0
	ovr=0.0
	vm=0.0
	xi=lxi
	yi=lyi
	zi=lzi
	rminnnsq=dn*dn
	rminnwsq=(dn/2.0+dw/2.0)*(dn/2.0+dw/2.0)
	! Interaction with positive ions
	do k=1,np
	rxi=xi-xp(k)
	ryi=yi-yp(k)
	rzi=zi-zp(k)
	rxi=rxi-L*anint(rxi/L)
	ryi=ryi-L*anint(ryi/L)
	rminpnsq=(dp(k)/2.0+dn/2.0)*(dp(k)/2.0+dn/2.0)
	dri2=rxi*rxi+ryi*ryi+rzi*rzi
	if(dri2>=rminpnsq) then
	dri=sqrt(dri2)
	vp=vp+(qn*qp(k)*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 600
	end if
	end do
	! Interaction with positively charged sheets
	do k=1,np
	zdis=abs(zi-zp(k))
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vps=vps+(qn*qp(k)*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	end do
	! Interaction with negative ions
	do k=1,nn
	if(k/=ind) then
	rxi=xi-xn(k)
	ryi=yi-yn(k)
	rzi=zi-zn(k)
	rxi=rxi-L*anint(rxi/L)
	ryi=ryi-L*anint(ryi/L)
	dri2=rxi*rxi+ryi*ryi+rzi*rzi
	if(dri2>=rminnnsq) then
	dri=sqrt(dri2)
	vn=vn+(qn*qn*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 600
	end if
	end if

	end do
	! Interaction with negatively charged sheets
	do k=1,nn
	zdis=abs(zi-zn(k))
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vns=vns+(qn*qn*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	end do
	! Interaction with the walls
	if(zi<(dn/2.0)) then
	ovr=1.0
	goto 600
	else
	if(zi>(W-dn/2.0)) then
	ovr=1.0
	goto 600
	else
	do k=1,nw
	rxi=xi-xw(k)
	ryi=yi-yw(k)
	rzi=zi-zw(k)
	rxi=rxi-L*anint(rxi/L)
	ryi=ryi-L*anint(ryi/L)
	dri2=rxi*rxi+ryi*ryi+rzi*rzi
	if(dri2>=rminnwsq) then
	dri=sqrt(dri2)
	vw=vw+(qn*qw*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 600
	end if
	end do
	end if
	end if
	! Interaction with wall charged sheets
	zdis=abs(zi)
	r1aux=0.5+(zdis/L)*(zdis/L)
	r1=sqrt(r1aux)
	r2aux=0.25+(zdis/L)*(zdis/L)
	r2=sqrt(r2aux)
	pos=(4.0*zdis*r1)/L
	arctan=atan(pos)
	EW=2.0*pi-4.0*arctan
	aux=log((0.5+r1)/r2)
	fi=-2.0*pi*zdis-4.0*L*aux+zdis*EW
	vws=nw*(qn*qw*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	vinti=vp+vps+vn+vns+vw+vws
	600 RETURN
	END SUBROUTINE energyneg
	! --------------------------------------------------------------------------
	SUBROUTINE writfinal(consist,vtot,vend,avv)

	USE mixdat
	IMPLICIT NONE
	real(kind=8),intent(in)::consist,vtot,vend,avv
	integer::j
	! Opening the file
	open(unit=8,file=fout2,status='unknown',form='formatted')
	! Saving the data
	write(unit=8,fmt=700)consist,vtot,vend,avv
	do j=1,npaux1
	write(unit=8,fmt=710)xp(j),yp(j),zp(j)
	end do
	write(unit=8,fmt=720)avv
	do j=npaux1+1,npaux2
	write(unit=8,fmt=710)xp(j),yp(j),zp(j)
	end do
	write(unit=8,fmt=720)avv
	do j=npaux1+2,np
	write(unit=8,fmt=710)xp(j),yp(j),zp(j)
	end do
	write(unit=8,fmt=720)avv
	do j=1,nn
	write(unit=8,fmt=710)xn(j),yn(j),zn(j)
	end do
	write(unit=8,fmt=720)avv
	do j=1,nw
	write(unit=8,fmt=710)xw(j),yw(j),zw(j)
	end do
	write(unit=8,fmt=720)avv
	do j=1,nhist
	write(unit=8,fmt=730)zcoord(j),cpav1(j),cpav2(j),cpav3(j),cnav(j)
	end do
	700 format(1x,E14.6,3x,E14.6,3x,E14.6,3x,E14.6)
	710 format(1x,E14.6,3x,E14.6,3x,E14.6)
	720 format(1x,E14.6)
	730 format(1x,E14.6,3x,E14.6,3x,E14.6,3x,E14.6,3x,E14.6)
	end file(unit=8)
	close(unit=8)
	RETURN
    END SUBROUTINE writfinal
	! -----------------------------------------------------------------------
	SUBROUTINE histogram(printnum)
	USE mixdat
	IMPLICIT NONE

	integer,intent(in)::printnum
	integer::j,k
	real::aux
	real::indexp(np),indexn(nn)
	real(kind=8)::concp1(nhist),concp2(nhist),concp3(nhist)
	real(kind=8)::concn(nhist)
	! Opening a file
	open(unit=7,file=fout3,status='unknown',form='formatted')
	! Initializing variables
	Do j=1,nhist
	concp1(j)=0.0
	concp2(j)=0.0
	concp3(j)=0.0
	concn(j)=0.0
	end do
	Do j=1,np
	indexp(j)=0.0
	End do
	Do j=1,nn
	indexn(j)=0.0
	End do
	Do k=1,nhist !Z coordinate
	Do j=1,npaux1 ! First kind of counterions
	If(indexp(j)==0.0) then
	aux=aint(zp(j)/deltaz)
	if(aux<=k) then
	indexp(j)=real(k)
	concp1(k)=concp1(k)+1.0
	end if
	end if
	end do
	Do j=1,npaux1 ! First kind of counterions
	If(indexp(j)==0.0) then
	aux=aint(zp(j)/deltaz)
	if(aux<=k) then
	indexp(j)=real(k)
	concp1(k)=concp1(k)+1.0
	end if
	end if
	end do
	Do j=npaux1+1,npaux2 ! Second kind of counterions
	If(indexp(j)==0.0) then
	aux=aint(zp(j)/deltaz)
	if(aux<=k) then
	indexp(j)=real(k)
	concp2(k)=concp2(k)+1.0
	end if
	end if
	end do
	Do j=npaux2+1,np ! Third kind of counterions
	If(indexp(j)==0.0) then
	aux=aint(zp(j)/deltaz)
	if(aux<=k) then
	indexp(j)=real(k)
	concp3(k)=concp3(k)+1.0
	end if

	end if
	end do
	Do j=1,nn ! Coions
	If(indexn(j)==0.0) then
	aux=aint(zp(j)/deltaz)
	if(aux<=k) then
	indexn(j)=real(k)
	concn(k)=concn(k)+1.0
	end if
	end if
	end do
	end do
	! Accumulating average concentrations
	Do k=1,nhist
	cpav1(k)=cpav1(k)+concp1(k)
	cpav2(k)=cpav2(k)+concp2(k)
	cpav3(k)=cpav3(k)+concp3(k)
	cnav(k)=cnav(k)+concn(k)
	End do
	! Saving the data
	If(mod(printnum,1000==0)) then
	    write(unit=7,fmt=800)printnum
	    Do k=1,nhist
		write(unit=7,fmt=810)zcoord(k),concp1(k), concp2(k), concp3(k), concn(k)
		  End do
	End if
	800 format(1x,I8)
	810 format(1x,E14.6,3x, E14.6,3x, E14.6,3x, E14.6,3x, E14.6)
	if(printnum==nstep) then
	end file (unit=7)
	close(unit=7)
	end if
	RETURN
	END SUBROUTINE histogram