program Vcpp ! Computes voltage due to alpha-helix ! of NR Rsat (x,y,0) and (x,y,-t) averaged over ! disk of radius 1 nm. Alpha-helix case. use defs; implicit none integer(i4b)::j,k,n,Nrandom,NR,ntries integer(i4b),parameter::NL=100 real(dp)::Vwsum,Vcsum,Vw,Vw0,Vwl,Vc,x,y,beg,spot real(dp)::rho,rhos,Vcbar,Vwbar real(dp),dimension(2)::rdn real(dp),parameter::elec=1.60217653D-19 ! C e's charge SI units real(dp),parameter::epsl=2.0_dp ! lipid dielectric real(dp),parameter::epsw=80.0_dp ! water dielectric (at 37 C?) real(dp),parameter::epsc=80.0_dp ! cytosol dielectric (at 37 C?) real(dp),parameter::eps0=8.854187817D-12 ! vacuum SI units real(dp),parameter::epswl = 0.5d0*(epsw+epsl),epsbar = epswl real(dp),parameter::epslc = 0.5d0*(epsl+epsc) real(dp),parameter::p=(epsw-epsl)/(epsw+epsl),pp=p,ppp=p*pp real(dp),parameter::t=5.0d-9 ! m write(6,*) 'Nrandom = ?' read(5,*) Nrandom call random_seed(); open(7,file='cppDV') ! Now compute the voltage at (x,y,0) due to ! the static R^Nsat (x,y,0) using Vw of (23=A49) ! and Vc of (23=A49). do NR = 5, 12 ntries = 0; Vwsum = 0.0d0; Vcsum = 0.0d0 do k = 1, Nrandom 10 call random_number(rdn); x= rdn(1); y = rdn(2) if (x**2+y**2 > 1.0d0) then go to 10 end if x = x*1.0d-9; y = y*1.0d-9; ntries = ntries + 1 do j = 1, NR beg = (0.5d0*0.16d-9)*real(NR,dp) spot = beg - (0.16d-9)*real(j,dp) rhos = x**2 + (y-spot)**2; rho = sqrt(rhos) Vw0 = (1.d0+p)/rho; Vwl = 0.0d0 do n = 1, NL Vwl = Vwl + (ppp**(n-1))/sqrt(rhos+(2*n*t)**2) end do Vw = Vw0 - (epsw*epsl/epswl**2)*pp*Vwl Vwsum = Vwsum + Vw Vc = 0.0d0 do n = 0, NL ! Vc of (23=A49) Vc = Vc + (ppp**n)/sqrt(rhos+((2*n+1)*t)**2) end do Vcsum = Vcsum + Vc end do end do Vwsum = Vwsum*elec/(4.0_dp*pi*eps0*epsw) ! V Vcsum = Vcsum*elec*epsl/(4.d0*pi*eps0*epswl*epslc)!V Vwbar = Vwsum/real(ntries,dp) Vcbar = Vcsum/real(ntries,dp) write(7,'(a,i3)')'NR =',NR write(7,'(a,f10.6,a)')' =',Vwbar,' V' write(7,'(a,f10.6,a)')' =',Vcbar,' V' write(7,'(a,f10.6,a)')'dVCPP =',Vcbar-Vwbar,' V' end do close(7) end program Vcpp