CODE ------- RESOLUTION OF THE POISSON-BOLTZMANN + PMF EQUATION WITH BOUNDARY CONDITIONS ON THE ELECTRIC FIELD --------------------- C by Carles Calero C Version: 05/01/2011 C C DESCRIPTION: This program solves the PB+PMF differential equation with boundary conditions y'(0) = -4*pi*0.71*DL*Sigma, y'(L) = 0 C (where y is the electrostatic potential). C - It uses the shooting method. That is, it actually solves the differential equation with a guess for the value of the potential at z=0 C and then checks whether the boundary condition at the upper boundary (z=L) is satisfied. C C - The program consists of three main blocks: C Block 1: Iterative process to find the potential at the origin that allows to solve the PB+PMF equation C (with boundary conditions on the electric field) C Block 2: Resolution of the PB+PMF equation using the "shooting" method C Block 3: Integration of the differential equation with boundary conditions at the lower boundary C C OBSERVATIONS: the PB+PMF differential equation is introduced in the subroutine "Derivades". C c INPUT: The user needs to introduce in the subroutine PB the parameters of the problem and integration: C - Debye length in nm C - Surface charge density in e/nm^2. C - Domain of integration C - Number of divisions of the domain of integration C When running the program, the user needs to introduce a guess for the value of the electrostatic potential at the origin. C C OUTPUT: "outfilePB+PMF.dat": five columns: position (nm), electrostatic potential (kT/e), electric field (kT/(e nm)), C concentration coions (M) and concentration counterions (M) C C C C C -- BLOCK 1: Iteration to find the potential at the origin that allows to solve the PB+PMF equation (with boundary conditions on the electric field)----- C C In this block, the guess for the value of the potential at the lower boundary (y0) is introduced by the user. In case y0 is close enough to the actual value of the potential at z=0 the differential equation will be solvable and the subroutine PB will adjust y0 until it finds the potential y0 which makes the solution to satisfy the boundary condition in bulk (y'(L)=0). In case y0 is not that close, the subroutine PB will not converge. It turns out, though, that in most cases, the closer y0 is to the correct value of y0 the further the integration goes. We make use of this observation to determine the increment dy that will provide a new guess for phi0 (done in the PB subroutine). C C IMPLICIT DOUBLE PRECISION (a-h, o-z) DIMENSION y(2, 0:1000000), yp(2) PARAMETER (pi= 3.14159265d0) COMMON DL, tfinold, dy, dy2 open(unit=11, file='main.log') C Introduce a guess for the potential at the origin. WRITE(*,*) 'Introduce guess for the potential at the origin:' READ(*,*) Y0 tfinold = 0.0 C If Y0 is not good enough, the program tries Y0+dy for the potential at the origin dy = 0.1 dy2 = 0.05 yp(1) = y0-dy yp(2) = yp(1)+dy2 num = 0 DO WHILE(num.LT.1000) num = num + 1 yp(1) = yp(1) + dy yp(2) = yp(1) + abs(dy2) WRITE(*,*) 'num=',num, ' yp1= ', yp(1), ' yp2= ', yp(2) write(*,*) ' dy=',dy, ' dy2=',dy2 WRITE(11,*) 'num=',num, ' yp1= ', yp(1), ' yp2= ', yp(2) CALL PB(yp(1), yp(2)) ENDDO WRITE(*,*) 'NUM= ', num close(11) END C ----- Block 2: Resolution of the PB+PMF equation using the "shooting" method ----------------------------- C C C Shooting method applied to the PB equation C Boundary conditions: y'(0) = -4*pi*0.71*DL*Sigma, y'(L) = 0 C SUBROUTINE PB(yp1, yp2) IMPLICIT DOUBLE PRECISION (a-h, o-z) DIMENSION y(2, 0:1000000), yp(2), fun(2) PARAMETER (pi= 3.14159265d0) COMMON DL, tfinold, dy, dy2 OPEN(UNIT=20,FILE='outfilePB+PMF.dat') WRITE(20,'(A,13X, A,13X, A,9X, A,8X, A,9X, A)')'#', 'z (nm)', 1 'V (kT/e)','E (kT/(e nm))', 'Ccoions (M)', 'Ccontra(M)' WRITE(20,*) ' ' OPEN(UNIT=21,FILE='PB2.log') C ..... Introduction of Parameters ............. !Introduce parameters of the problem: DL= 25.5 ! Debye length (in nm) Sigma = 0.1 ! surface charge density (in e/nm²) !Introduce parameters of the integration: dintegration = 6.0*DL ! length of the domain of integration ninterv = 10000 ! Number of divisions of the interval h= dintegration/ninterv ! integration step t0=0 ! origin neqs=2 ndim=2 y(2,0)= - 4*pi*0.71*Sigma error=1.d-4 tfin = 0.0 C ............................................ yp(1) = yp1 yp(2) = yp2 C The values for y and y' at the lower boundary are introduced DO icas=1,2 y(1,0) = yp(icas) c Integration CALL EQDIFS(t0,h,ninterv,neqs,ndim,y,error, tfin) WRITE(*,*) 'Value at origin =', Y(1,0), 'YB=', Y(2,ninterv) fun(icas) = y(2,ninterv)- 0.0 C Determination of the values of the increments dy and dy2 IF(icas.EQ.1) THEN IF((tfin.GT.0.0).AND.(tfin.GT.tfinold)) THEN dy = dy tfinold = tfin EXIT ELSE IF((tfin.GT.0.0).AND.(tfin.LT.tfinold)) THEN dy = -dy/2.0 dy2 = -dy2/2.0 tfinold = tfin EXIT ENDIF IF(isnan(y(2,ninterv))) EXIT ELSE IF(icas.EQ.2) THEN IF(tfin.GT.0.0) THEN dy2 = dy2/2.0 dy = 0.0 EXIT ENDIF ENDIF ENDDO IF(tfin.NE.0.0) GOTO 200 c Interpolation: Regula Falsi 100 ypnova=(yp(1)*fun(2)-yp(2)*fun(1))/(fun(2)-fun(1)) y(1,0) = ypnova CALL EQDIFS(t0,h,ninterv,neqs,ndim,y,error,tfin) IF(tfin.NE.0.0) GOTO 200 WRITE(*,*) 'Valor origen=', Y(1,0), 'YB=', Y(2,ninterv) WRITE(21,*) 'Valor origen=', Y(1,0), 'YB=', Y(2,ninterv) funnova = y(2,ninterv)-0.0 write(*,*) 'funnova= ', funnova write(21,*) 'funnova= ', funnova write(*,*) 'D2' if(isnan(funnova)) goto 200 write(*,*) 'D3' IF(ABS(funnova).LT.1.D-5) THEN icont = 0 ELSE icont=1 ENDIF IF (icont.EQ.0) THEN DO I=0,ninterv WRITE(*,'(1x,3F20.8)') real(I)*h, y(1,i), y(2,i) WRITE(20,'(1x,5F20.8)') real(I)*h, y(1,i), y(2,i), 1 (1.0/(8.0*3.1416*0.6022*DL*DL))*EXP(-y(1,i)), 1 (1.0/(8.0*3.1416*0.6022*DL*DL))*EXP(y(1,i)-1.0* 1 (-real(I)*h*real(I)*h*EXP(-13.1621*(real(I)*h-0.614988)))) ENDDO STOP ENDIF c New values IF (ABS(yp(1)-ypnova) .LT. ABS(yp(2)-ypnova)) THEN yp(2)= ypnova fun(2) =funnova ELSE yp(1) = ypnova fun(1) = funnova ENDIF GOTO 100 200 CLOSE(20) CLOSE(21) RETURN END C ----- Block 3: Integration of the differential equation with boundary conditions at the lower boundary ----------------------------- SUBROUTINE EQDIFS(t0,h,np,neqs,ndim,y,error,tfin) C C C C C C IMPLICIT DOUBLE PRECISION (a-h, o-z) PARAMETER (nemax=10) DIMENSION y(ndim, 0:np) DIMENSION y1(nemax), yn(nemax) IF(neqs.GT.nemax) THEN WRITE(*,*) 'Too many equations in EQDIFS' STOP ENDIF hh = h NVOLTES =1 t = t0 DO I = 1, NP 100 ierr = 0 t1 = t DO 10 L = 1, neqs 10 y1(L) = y(L,I-1) DO K=1,NVOLTES CALL RK(t1,hh,neqs,y1,yn,error,ierror) C Step failed IF (ierror.EQ.2) THEN hh=hh/2.0 NVOLTES = 2*nvoltes IF (nvoltes.GT.10000) THEN write(*,*) 'A2' tfin = t1 c WRITE(*,*) 'Too many subdivisions at EQDIFS', '; t = ', t1, EXIT ENDIF GOTO 100 ENDIF IF(tfin.NE.0.0) EXIT C Correct step ierr = MAX(ierr,ierror) t1 = t1+hh DO 20 L=1, neqs 20 y1(L) = yn(L) ENDDO IF(tfin.NE.0.0) EXIT t = t+h DO 30 L=1,neqs 30 y(L,I)=yn(L) IF(ierr.EQ.0.AND.nvoltes.GT.1) THEN hh=hh*2 nvoltes=nvoltes/2 ENDIF ENDDO RETURN END CODI ---------------------- EULER----------------------------- SUBROUTINE EULER 1 (t, h, neqs, y, ynou, error, ierror) C C C IMPLICIT DOUBLE PRECISION (a-h, o-z) PARAMETER (nemax=10, nsub=6) DIMENSION Y(neqs), YNOU(neqs) DIMENSION np(0:nsub), YY(0:nsub,nemax), hh(0:nsub) DIMENSION f(nemax,0:16), yerror(nemax), der(nemax), 1 der0(nemax) c EXTERNAL derivades data np/1,2,3,5,7,10,16/ IF (neqs.GT.nemax) THEN WRITE(*,*) "Too many equations in EULER " STOP ENDIF C DEPARTING POINT IERROR=0 DO I = 1, neqs f(I,0) = Y(I) ENDDO C STEP H hh(0) = h CALL derivades(t, y, der0) DO I = 1, neqs yy(0,I) = y(I) + h*der0(I) ENDDO C Calculations for Richardson's extrapolation DO 20 J=1, nsub npunts = np(J) h1 = h/npunts hh(J) = h1 C The first point (0) is already known DO I =1, neqs f(I,1)= f(I,0) + h1*der0(I) ENDDO DO 10 I =2,npunts x = t + (i-1)*h1 CALL derivades(x,f(1,i-1),der) DO K=1, neqs f(K,I) = f(K,I-1) + h1*der(K) ENDDO 10 CONTINUE DO I =1, neqs yy(J,I) = f(I, npunts) ENDDO C Extrapolation xerror = 0 DO I = 1,neqs CALL interpolacio 1 (hh,yy(0,I), j+1,0.d0, ynou(I), yerror(I)) xerror=MAX(xerror,ABS(yerror(I)/ynou(I))) ENDDO C Test on the error and assessment of the covergence IF (xerror.LE.error) THEN IF(J.LE.3) ierror =0 IF(J.GT.3) ierror =1 RETURN ENDIF 20 CONTINUE C ENDS WITH NO CONVERGENCE ierror = 2 RETURN END CODI ------------------- SUBROUTINE RUNGE-KUTTA ---------------------------------------- SUBROUTINE RK 1 (t,h,neqs,y,ynou,error,ierror) c Integration using the Runge-Kutta Cash-Karp method c W.H.Press, Comp. Phys. 6 (1992) 188 c t -> Origin c h -> Integration step c NEQS -> Number of equations c Y(1..NEQS) -> Initial values of the function c YNOU -> Final value of the function c ERROR -> Relative tolerance error c IERROR -> 0 Convergence, short step c -> 1 Convergence c -> 2 It doesn't converge (very long integration step) c DERIVADES -> Subroutine EXTERNAL.Calculates DERIVADES IMPLICIT REAL*8 (a-h,o-z) PARAMETER (nemax=10) PARAMETER 1 (a2=0.2d0,a3=0.3d0,a4=0.6d0,a5=1.d0,a6=0.875d0) PARAMETER(b21=0.2d0,b31=3.d0/40.d0,b41=0.3d0, 1 b51=-11.d0/54.d0,b61=1631.d0/55296.d0) PARAMETER(b32=9.d0/40.d0,b42=-0.9d0,b52=2.5d0, 1 b62=175.d0/512.d0) PARAMETER 1 (b43=1.2d0,b53=-70.d0/27.d0,b63=575.d0/13824.d0) PARAMETER(b54=35.d0/27.d0,b64=44275.d0/110592.d0) PARAMETER(b65=253.d0/4096.d0) PARAMETER(c1=37.d0/378.d0,c3=250.d0/621.d0, 1 c4=125.d0/594.d0, c6=512.d0/1771.d0) PARAMETER(d1=2825.d0/27648.d0,d3=18575.d0/48384.d0, 1 d4=13525.d0/55296.d0,d5=277.d0/14336.d0,d6=0.25d0) DIMENSION y(neqs),ynou(neqs) DIMENSION xk1(nemax),xk2(nemax),xk3(nemax), 1 xk4(nemax),xk5(nemax),xk6(nemax), yy(nemax) c CALL derivades(t,y,xk1) DO 10 l=1,neqs 10 yy(l)=y(l)+h*b21*xk1(l) CALL derivades(t+a2*h,yy,xk2) DO 20 l=1,neqs 20 yy(l)=y(l)+h*(b31*xk1(l)+b32*xk2(l)) CALL derivades(t+a3*h,yy,xk3) DO 30 l=1,neqs 30 yy(l)=y(l)+h*(b41*xk1(l)+b42*xk2(l)+b43*xk3(l)) CALL derivades(t+a4*h,yy,xk4) DO 40 l=1,neqs 40 yy(l)=y(l)+ 1 h*(b51*xk1(l)+b52*xk2(l)+b53*xk3(l)+b54*xk4(l)) CALL derivades(t+a5*h,yy,xk5) DO 50 l=1,neqs 50 yy(l)=y(l)+ 1 h*(b61*xk1(l)+b62*xk2(l)+b63*xk3(l)+b64*xk4(l)+ 1 b65*xk5(l)) CALL derivades(t+a6*h,yy,xk6) err=0 DO 60 l=1,neqs ynou(l)=y(l)+ 1 h*(c1*xk1(l)+c3*xk3(l)+c4*xk4(l)+c6*xk6(l)) c Error yy(l)=h*((c1-d1)*xk1(l)+(c3-d3)*xk3(l)+ 1 (c4-d4)*xk4(l)- d5*xk5(l)+(c6-d6)*xk6(l)) err= 1 MAX(error,2*ABS(yy(l))/(ABS(y(l))+ABS(ynou(l)))) 60 CONTINUE IF (err .LE. error) THEN ierror=1 IF (error/err .GT. 10.d0) ierror=0 ELSE ierror=2 ENDIF RETURN END CODI ------------------------- SUBROUTINE derivades ------------------------- SUBROUTINE derivades(t, y, der) implicit DOUBLE PRECISION (a-h, o-z) COMMON DL PARAMETER (pi= 3.14159265d0) dimension y(2), der(2) c der(1) := y' der(1)= y(2) c der(2):= y'' der(2)=-(1.0/(2.0*DL*DL))*(EXP(-y(1))-EXP(y(1)-1.0* 1 (-t*t*EXP(-13.1621*(t-0.614988))))) return end CODI --------------------- SUBROUTINE INTERPOLACIO -------------------------- SUBROUTINE INTERPOLACIO(xa,ya,N,X,Y,ERROR) implicit DOUBLE PRECISION (a-h, o-z) PARAMETER (nmax = 20) DIMENSION xa(n), ya(n), d(nmax,2) IF(n.GT.nmax) THEN WRITE(*,*) 'Error subroutine INTERPOLACIO, n>', nmax STOP ENDIF c Starts iteration DO 10 I=1,n 10 d(I,1) = ya(I) c Iteration in = 1 io = 2 DO 20 j=1,n-1 DO 30 i=1,n-j den=xa(i+j) -xa(i) IF(den.EQ.0.d0) THEN WRITE(*,*) 'Error interpolation', 1 ' two identical abscisses' STOP ENDIF 30 d(i,io) = ((x-xa(i))*d(i+1,in)-(x-xa(i+j))*d(i,in))/den ix=in in=io 20 io=ix C Result y=d(1,in) C Estimate of error error = MIN(ABS(y-d(1,io)),ABS(y-d(2,io))) RETURN END logical function isnan(x) double precision x isnan = x .ne. x return end