mentecaro.for

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

										 MODULE mixdat
	! Declaration of global parameters and global variables
	IMPLICIT NONE
	SAVE
	integer,parameter::nstep=1.8E5,nn=324,nw=36
	integer,parameter::npaux1=40,npaux2=100,np=220
	real,parameter::eps=78.5,epso=8.8541E-22,elec=1.60219E-19,pi=3.1415297
	real,parameter::L=30.0,W=600.0,dn=4.25,dw=3.0
	real,parameter::deltaz=5.0
	real,parameter::kb=1.38062E-23,temp=298
	real,parameter::qn=-1.0,qw=-1.0
	integer::nhist !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	real,dimension(np)::dp,qp
	real,dimension(np)::xp,yp,zp
	real,dimension(nn)::xn,yn,zn
	real,dimension(nw)::xw,yw,zw
	real, dimension(nhist)::zcoord,cpav1,cpav2,cpav3,cnav
	character (len=256):: fout1, fout2, fout3
	END MODULE mixdat
	! -----------------------------------------------------------------------
	PROGRAM Predlmix
	! Monte Carlo simulation of the EDL formation in mixtures of electrolytes
	! NVT Ensemble
	! Patricia Taboada Serrano (2004)
	USE mixdat
	IMPLICIT NONE
	! Declaration of variables
	integer::step,correc,iprint,icalc,limit,i,j,ihisto,isave
	integer::icount,fcount,index
	real::ovrlap,ovrold,ovrnew
	real::countc,acm,acatma,ratio,cont
	real::drmax,dummy,dummy1,dummy2,dummy3
	real::xiold,yiold,ziold,xinew,yinew,zinew
	real::rminppsq,rminnnsq,rminpnsq,rminpwsq,rminnwsq
	real::rminpmsq,rminnmsq,rminwmsq
	real(kind=8)::avv,acv,consist,vend,vtot
	real(kind=8)::aux,aux1,vold,vnew,deltv
	real:: p1,p2
	namelist /forcein / p1,p2, fout1, fout2, fout3
	read(5,forcein)
	! Initialization of variables
	drmax=1.5*dw !(ion displacement)
	correc=10 !(correction of displacement)
	iprint=1000 !(printout step)

	icalc=10 !(calculation of force)
	limit=1.7E5 !(equilibration period)
	CALL random_seed() !Initializing random number generator
	do i=1,npaux1
	dp(i)=9.0
	qp(i)=3.0
	end do
	do i=npaux1+1,npaux2
	dp(i)=6.0
	qp(i)=2.0
	end do
	do i=npaux2+1,np
	dp(i)=4.25
	qp(i)=1.0
	end do
	do i=1,np
	xp(i)=0.0
	yp(i)=0.0
	zp(i)=0.0
	end do
	do i=1,nn
	xn(i)=0.0
	yn(i)=0.0
	zn(i)=0.0
	end do
	! Initialization of Z coordinate and conc. Profiles
	nhist=int(W/deltaz)
	zcoord(1)=deltaz
	Do j=2,nhist
	zcoord(j)=zcoord(j-1)+deltaz
	end do
	Do j=1,nhist
	cpav1(j)=0.0
	cpav2(j)=0.0
	cpav3(j)=0.0
	cnav(j)=0.0
	end do
	! Equispaced ions uniformly distributed on the surface
	! Ions being numbered from left upper corner
	do i=1,nw
	zw(i)=0.0
	end do
	j=1
	do i=1,nw
	xw(i)=-12.5+(j-1)*(dw+(L-6.0*dw)/6.0)
	j=j+1
	if(j>6) then
	j=1
	end if
	end do
	icount=1
	fcount=6
	do j=1,6
	if(j>1) then
	icount=6*j-5
	fcount=icount+6

	end if
	do i=icount,fcount
	yw(i)=12.5-(j-1)*(dw+(L-6.0*dw)/6.0)
	end do
	end do
	! Generation of Initial Configuration
	CALL iniconf(aux)
	CALL writini(aux)
	! Initialization of accumulators
	vtot=0.0 !(total potential energy)
	acv=0.0 !(acummulated average potential energy)
	vend=0.0 !(final potential energy)
	acm=0.0 !(total number of attempted moves)
	acatma=0.0 !(total number of accepted moves)
	cont=0.0 !(# of cycles over which average is taken)
	avv=0.0 !(average potential energy)
	ovrlap=0.0 !(overlap criterium)
	ratio=0.0 !(# accepted moves/# total moves)
	countc=0.0 !(# steps over which conc. profiles are averaged)
	! Calculation of initial energy
	CALL sumup(ovrlap,vtot)
	write(*,fmt=110)vtot,ovrlap
	110 format(1x,E14.6,3x,f10.6)
	if(ovrlap==1.0) then
	goto 100
	end if
	! Monte Carlo simulation
	DO step=1,nstep ! Loop over all cycles
	! Movement of positive ions
	DO i=1,np
	index=i
	xiold=xp(i)
	yiold=yp(i)
	ziold=zp(i)
	vold=0.0
	vnew=0.0
	ovrold=0.0
	ovrnew=0.0
	! Energy of ion i in the old configuration
	CALL energyp(xiold,yiold,ziold,index,vold,ovrold)
	! Move i and generate image ions
	CALL random_number(dummy1)
	xinew=xiold+(2.0*dummy1-1.0)*drmax
	CALL random_number(dummy2)
	yinew=yiold+(2.0*dummy2-1.0)*drmax
	CALL random_number(dummy3)
	zinew=ziold+(2.0*dummy3-1.0)*drmax
	xinew=xinew-L*anint(xinew/L)
	yinew=yinew-L*anint(yinew/L)
	! Energy of ion i in the new configuration
	CALL energyp(xinew,yinew,zinew,index,vnew,ovrnew)
	! Check movement acceptance
	deltv=(vnew-vold)/(kb*temp)
	if(ovrnew==0.0) then
	if(deltv<0.0) then
	vtot=vtot+deltv*kb*temp
	xp(i)=xinew
	yp(i)=yinew
	zp(i)=zinew
	acatma=acatma+1.0
	else
	CALL random_number(dummy)
	if((exp(-deltv))>dummy) then
	vtot=vtot+deltv*kb*temp
	xp(i)=xinew
	yp(i)=yinew
	zp(i)=zinew
	acatma=acatma+1.0
	end if
	end if
	end if
	acm=acm+1.0
	ratio=acatma/acm
	! Accumulation of average energy
	if(step>limit) then
	acv=acv+vtot
	cont=cont+1.0
	avv=acv/cont
	end if
	end do ! End cycle over all positive ions
	! Movement of negative ions
	do i=1,nn
	index=i
	xiold=xn(i)
	yiold=yn(i)
	ziold=zn(i)
	vold=0.0
	vnew=0.0
	ovrold=0.0
	ovrnew=0.0
	! Energy of ion i in the old configuration
	CALL energyneg(xiold,yiold,ziold,index,vold,ovrold)
	! Move ion i and generate image ions
	CALL random_number(dummy1)
	xinew=xiold+(2.0*dummy1-1.0)*drmax
	CALL random_number(dummy2)
	yinew=yiold+(2.0*dummy2-1.0)*drmax
	CALL random_number(dummy3)
	zinew=ziold+(2.0*dummy3-1.0)*drmax
	xinew=xinew-L*anint(xinew/L)
	yinew=yinew-L*anint(yinew/L)
	! Energy of ion i in the new configuration
	CALL energyneg(xinew,yinew,zinew,index,vnew,ovrnew)

	! Checking for movement acceptance
	deltv=(vnew-vold)/(kb*temp)
	if(ovrnew==0.0) then
	if(deltv<0.0) then
	vtot=vtot+deltv*kb*temp
	xn(i)=xinew
	yn(i)=yinew
	zn(i)=zinew
	acatma=acatma+1.0
	else
	CALL random_number(dummy)
	if((exp(-deltv))>dummy) then
	vtot=vtot+deltv*kb*temp
	xn(i)=xinew
	yn(i)=yinew
	zn(i)=zinew
	acatma=acatma+1.0
	end if
	end if
	end if
	acm=acm+1.0
	ratio=acatma/acm
	! Accumulation of average energy
	if(step>limit) then
	acv=acv+vtot
	cont=cont+1.0
	avv=acv/cont
	end if
	end do ! End of cycle over all negative ions
	! Periodic operations
	if(mod(step,iprint)==0) then
	write(*,fmt=150)vtot,acatma,acm,ratio
	end if
	150 format(1x,E14.6,3x,E14.6,3x,E14.6,3x,E14.6)
	if(mod(step,correc)==0) then
	if(ratio>0.3) then
	drmax=drmax*1.05
	else
	drmax=drmax*0.95
	end if
	end if
	! Calculating force
	if(step>limit) then
	if(mod(step,icalc)==0) then
	CALL histogram(step)
	countc=countc+1.0
	end if
	end if
	END DO ! End of cycles
	! Average concentration profiles
	DO j=1,nhist
	cpav1(j)=cpav1(j)/countc
	cpav2(j)=cpav2(j)/countc
	cpav3(j)=cpav3(j)/countc

	cnav(j)=cnav(j)/countc
	END DO
	! Checking consistency of results
	CALL sumup(ovrlap,vend)
	consist=abs(vtot-vend)
	! Writing final and average values
	CALL writfinal(consist,vtot,vend,avv)
	100 STOP
	END PROGRAM
	! --------------------------------------------------------------------------
	SUBROUTINE iniconf(overlap)
	USE mixdat
	IMPLICIT NONE
	real,intent(out)::overlap
	integer::j,k
	real::drx,dry,drz,pprsq,pnrsq,pwrsq
	real::nnrsq,nwrsq
	real::dummy1,dummy2,dummy3
	real::rminppsq,rminnnsq,rminpnsq,rminpwsq,rminnwsq
	real  rminnmsq
	! Initialization of variables
	overlap=0.0
	rminnnsq=dn*dn
	rminnwsq=(dn/2.0+dw/2.0)*(dn/2.0+dw/2.0)
	! rminnmsq=(dn/2.0+dw/2.0)*(dn/2.0+dw/2.0)
	 rminnmsq=rminnwsq
	! Placing positive ions in the box
	do j=1,np
	200 overlap=0.0
	rminpwsq=(dp(j)/2.0+dw/2.0)*(dp(j)/2.0+dw/2.0)
	CALL random_number(dummy1)
	xp(j)=(0.5-dummy1)*L
	CALL random_number(dummy2)
	yp(j)=(0.5-dummy2)*L
	CALL random_number(dummy3)
	zp(j)=dummy3*(W-dp(j)/2.0)
	! Checking overlap with previously placed positive ions
	if(j>1) then
	do k=1,j-1
	drx=xp(j)-xp(k)
	dry=yp(j)-yp(k)
	drz=zp(j)-zp(k)
	drx=drx-L*anint(drx/L)
	dry=dry-L*anint(dry/L)
	rminppsq=(dp(j)/2.0+dp(k)/2.0)*(dp(j)/2.0+dp(k)/2.0)
	pprsq=drx*drx+dry*dry+drz*drz
	if(pprsq<rminppsq) then
	overlap=1.0
	goto 200
	end if
	end do
	end if

	! Checking overlap with the walls
	if(zp(j)<(dp(j)/2.0)) then
	overlap=1.0
	goto 200
	else
	if(zp(j)>(W-dp(j)/2.0)) then
	overlap=1.0
	goto 200
	else
	do k=1,nw
	drx=xp(j)-xw(k)
	dry=yp(j)-yw(k)
	drz=zp(j)-zw(k)
	drx=drx-L*anint(drx/L)
	dry=dry-L*anint(dry/L)
	pwrsq=drx*drx+dry*dry+drz*drz
	if(pwrsq<rminpwsq) then
	overlap=1.0
	goto 200
	end if
	end do
	end if
	end if
	! Placing negative ions in the box
	do j=1,nn
	210   overlap=0.0

	CALL random_number(dummy1)
	xn(j)=(0.5-dummy1)*L
	CALL random_number(dummy2)
	yn(j)=(0.5-dummy2)*L
	CALL random_number(dummy3)
	zn(j)=dummy3*(W-dn/2.0)
	! Checking overlap with previously placed negative ions
	if(j>1) then
	do k=1,j-1
	drx=xn(j)-xn(k)
	dry=yn(j)-yn(k)
	drz=zn(j)-zn(k)
	drx=drx-L*anint(drx/L)
	dry=dry-L*anint(dry/L)
	nnrsq=drx*drx+dry*dry+drz*drz
	if(nnrsq<rminnnsq) then
	overlap=1.0
	goto 210
	end if
	end do
	end if
	! Checking overlap with positive ions
	do k=1,np
	drx=xn(j)-xp(k)
	dry=yn(j)-yp(k)
	drz=zn(j)-zp(k)
	drx=drx-L*anint(drx/L)
	dry=dry-L*anint(dry/L)
	rminpnsq=(dp(k)/2.0+dn/2.0)*(dp(k)/2.0+dn/2.0)
	pnrsq=drx*drx+dry*dry+drz*drz
	if(pnrsq<rminpnsq) then
	overlap=1.0
	goto 210
	end if
	end do

	! Checking overlap with the walls
	if(zn(j)<(dn/2.0)) then
	overlap=1.0
	goto 210
	else
	if(zn(j)>(W-dn/2.0)) then
	overlap=1.0
	goto 210
	else
	do k=1,nw
	drx=xn(j)-xw(k)
	dry=yn(j)-yw(k)
	drz=zn(j)-zw(k)
	drx=drx-L*anint(drx/L)
	dry=dry-L*anint(dry/L)
	nwrsq=drx*drx+dry*dry+drz*drz
	if(nwrsq<rminnwsq) then
	overlap=1.0
	goto 210
	end if
	end do
	end if
	end if
	end do ! End placing the negative ions
	RETURN
	!end do   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	END SUBROUTINE iniconf
	! ---------------------------------------------------------------------------
	SUBROUTINE writini(ovr)
	USE mixdat
	IMPLICIT NONE
	real,intent(in)::ovr
	integer::k
	! Opening the file
	open(unit=9,file=fout1,status='unknown',form='formatted')
	! Saving the data
	write(unit=9,fmt=310)ovr
	do k=1,npaux1
	write(unit=9,fmt=320)xp(k),yp(k),zp(k)
	end do
	write(unit=9,fmt=310)ovr
	do k=npaux1+1,npaux2
	write(unit=9,fmt=320)xp(k),yp(k),zp(k)
	end do
	write(unit=9,fmt=310)ovr
	do k=npaux2+1,np
	write(unit=9,fmt=320)xp(k),yp(k),zp(k)
	end do
	write(unit=9,fmt=310)ovr

	do k=1,nn
	write(unit=9,fmt=320)xn(k),yn(k),zn(k)
	end do
	write(unit=9,fmt=310)ovr
	do k=1,nw
	write(unit=9,fmt=320)xw(k),yw(k),zw(k)
	end do
	310 format(1x,f10.6)
	320 format(1x,E14.6,3x,E14.6,3x,E14.6)
	end file(unit=9)
	close(unit=9)
	RETURN
	END SUBROUTINE writini
	! --------------------------------------------------------------------------------------------------
	SUBROUTINE sumup(ovr,vint)
	USE mixdat
	IMPLICIT NONE
	real(kind=8),intent(out)::ovr
	real,intent(out)::vint
	integer::i,j,k
	real(kind=8)::vpp,vpps,vpn,vpns,vpw,vpws
	real(kind=8)::vnn,vnps,vnns,vnw,vnws
	real::rxi,ryi,rzi,dri2,dri
	real::xi,yi,zi,zdis
	real::r1,r1aux,r2,r2aux,pos,arctan
	real::EW,aux,fi
	real::qpi,dpi
	real::rminppsq,rminnnsq,rminpnsq,rminpwsq,rminnwsq
	! Initialization of variables
	ovr=0.0
	vint=0.0
	vpp=0.0
	vpps=0.0
	vpn=0.0
	vpns=0.0
	vpw=0.0
	vpws=0.0
	vnn=0.0
	vnps=0.0
	vnns=0.0
	vnw=0.0
	vnws=0.0
	!vwsm=0.0
	rminnnsq=dn*dn
	rminnwsq=(dn/2.0+dw/2.0)*(dn/2.0+dw/2.0)
	! Sum over all positive ions
	do i=1,np
	xi=xp(i)
	yi=yp(i)
	zi=zp(i)
	dpi=dp(i)

	qpi=qp(i)
	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 of positive ions with positive ions
	if(i<np) then
	do k=i+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)
	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)
	vpp=vpp+(qpi*qp(k)*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 400
	end if
	end do
	end if
    ! Interaction of positive 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
	vpps=vpps+(qpi*qp(k)*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	end do
	! Interaction of positive ions 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)
	vpn=vpn+(qpi*qn*elec*elec)/(4.0*pi*eps*epso*dri)
	else
	ovr=1.0
	goto 400
	end if
	end do
	! Interaction of positive 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
	vpns=vpns+(qpi*qn*elec*elec*fi)/(4.0*pi*eps*epso*L*L)
	end do
	! Interaction of positive ions with the walls
	if(zi<(dpi/2.0)) then
	ovr=1.0
	goto 400