c ******************************************************************** c Brownian dynamics algorithm. taking into account the hydrodynamics. c Second-order, Ermak-McCammon c hydro0 -- modified scheme (updating D here) c hyd0 -- modifed scheme (use old D) c ******************************************************************** c brownian dynamics algorithm. taking into account the hydrodynamics. c revised second-order BD (while updating diffu_tensor) subroutine bd_hydro(n,ntime) implicit double precision(a-h, o-z) parameter(maxa=500,maxn=maxa*3,bolz=8.31d-3/4.18,one3=1.d0/3.d0) integer nout, nn,ntime,n3,n double precision gc(maxn),xold(maxa,3),gcnew(maxn),fcnew, + y(maxn),rran(maxn),fold(maxa,3),vold(maxa,3),uold(maxa,3) + , winv(maxn),ftor(maxa),phiold(maxa),sum(maxn) common/matrix/x(maxa,3),bangle(maxa) common/matrix2/ r(maxa,maxa) common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma, +amass,am(maxn),coef,dtime2 common/diffu_tensor/diffu(maxn,maxn),w(maxn,maxn),varr(maxn) common/myrandom/amean,var(maxn),fvar(maxn) common/euler/alpha(maxa),beta(maxa),gama(maxa),phi(maxa) common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3) common/factor/dc1,dc2,bead2r,beadr2r,beadf1,beadf2 common/ischeme/icircular,iele,ihydro,iclose,iknot,neuler c debug c write(*,*) 'dc1=',dc1,'dc2=',dc2 c write(*,*) 'bead2r=',bead2r,' beadr2r=',beadr2r c write(*,*) 'beadf1=',beadf1,'beadf2=',beadf2 c update r(i,j) n3=n*3 call update_rij(n) call force4(n,gc,ntime) do i = 1, n xold(i,1) = x(i,1) xold(i,2) = x(i,2) xold(i,3) = x(i,3) if (icircular .eq. 1) then do j = 1, 3 fold(i,j) = f(i,j) uold(i,j) = u(i,j) vold(i,j) = v(i,j) enddo phiold(i) = phi(i) endif enddo c update D c update dphiz(rz),f,v by ftor if ((icircular .eq. 1).and.(iclose.eq.1)) then iout = 0 call torsion_bd(n,ftor,phiold,iout) endif C initialize the upper half 3x3 matrix to be zero. do idi = 3, n3, 3 diffu(idi-2,idi-1) = 0.d0 diffu(idi-2,idi) = 0.d0 diffu(idi-1,idi) = 0.d0 enddo C mean value and variance for random variable y. C the diagonal element of diffu set to be KT/(6*pi*eta*a) c this part is only done once in circle.f do i = 1, n3 diffu(i,i) = dc1 enddo C call to get normalized random variable y. C =0 =2*dtime (when i=j) call rannv2(n3,y,amean,varr) do i = 1, n-1 idi = 3*i-3 do j = i+1, n idj = 3*j-3 rij = r(i,j) vrij = 1.d0/rij vr2 = vrij*vrij if (rij .gt. bead2r) then C non-overlapping beads do i1 = 1, 3 do j1 = i1, 3 if (j1 .eq. i1) then tmp = (x(i,i1)-x(j,i1))**2*vr2 tmp2 = beadr2r*vr2 diffu(idi+i1,idj+j1)=dc2*vrij*(1.d0+tmp2*one3+ + tmp*(1.d0-tmp2)) else tmp=(x(i,i1)-x(j,i1))*(x(i,j1)-x(j,j1))*vr2*vrij diffu(idi+i1,idj+j1)=dc2*tmp*(1.d0-beadr2r*vr2) diffu(idi+j1,idj+i1)=diffu(idi+i1,idj+j1) endif enddo enddo else C overlapping beads do i1 = 1, 3 do j1 = i1, 3 if (j1 .eq. i1) then tmp = x(i,i1)-x(j,i1) tmp = tmp*tmp*vrij diffu(idi+i1,idj+j1)=dc1*(1.d0-beadf1*rij+beadf2*tmp) else tmp=(x(i,i1)-x(j,i1))*(x(i,j1)-x(j,j1))*vrij diffu(idi+i1,idj+j1)=dc1*tmp*beadf2 diffu(idi+j1,idj+i1)=diffu(idi+i1,idj+j1) endif enddo enddo endif enddo enddo C assemble the random force using the procedure C described in Ermak-McCammon's paper. do i = 1, n3 rran(i) = 0.d0 do j = 1, i summ = 0.d0 do k = 1, j-1 summ = summ + w(i,k)*w(j,k) enddo tmp10 = diffu(j,i)-summ if (j .eq. i) then if (tmp10 .lt. 0.d0) then print *, "diffu negative", diffu(j,i)-summ call bd_stop(n) stop endif w(i,j)=dsqrt(tmp10) winv(j) = 1.d0/w(i,j) else w(i,j)=tmp10*winv(j) c w(i,j)=tmp10/w(j,j) endif rran(i) = rran(i) + w(i,j)*y(j) enddo enddo C Generate the new position vector. do i = 1, n3 sum(i) = 0.d0 enddo do i=1, n3 do j = 1, n3 if (i .ge. j) then sum(i) = sum(i) + diffu(j,i)*gc(j) else sum(i) = sum(i) + diffu(i,j)*gc(j) endif enddo enddo c in my model i=3*id+jd-3 jd=1:3 id=1:n ii=1 do i=1,n do j=1,3 x(i,j)=x(i,j)+dtime*sum(ii)+rran(ii) ii=ii+1 enddo enddo call update_rij(n) call update_vector(n) if ((icircular .eq. 1).and.(iclose.eq.1)) then call update_euler(n) endif C Go further to use second-order BD. call force4(n,gcnew,ntime) if ((icircular .eq. 1).and.(iclose.eq.1)) then do i = 1, n do j = 1, 3 f(i,j) = fold(i,j) v(i,j) = vold(i,j) u(i,j) = uold(i,j) enddo enddo iout = 1 call torsion_bd(n,ftor,phiold,iout) endif do i = 1, n3 sum(i) = 0.d0 enddo do i=1, n3 do j = 1, n3 if (i .ge. j) then sum(i) = sum(i) +diffu(j,i)*(gc(j)+gcnew(j)) else sum(i) = sum(i) +diffu(i,j)*(gc(j)+gcnew(j)) endif enddo enddo c in my model i=3*id+jd-3 jd=1:3 id=1:n ii=1 do i=1,n do j=1,3 x(i,j)=xold(i,j)+0.5*dtime*sum(ii)+rran(ii) ii=ii+1 enddo enddo call update_rij(n) call update_vector(n) if ((icircular .eq. 1).and.(iclose.eq.1)) then call update_euler(n) endif c write(990,*) 0.5d0*dtime*sum,coef*gc(n3),ntime c write(991,*) rran(n3),y(n3),ntime 777 format(g9.5,1x,g9.5,1x,g9.5,1x,g9.5,1x,i4) c do i = 1, n3 c do j = 1, n3 c if (i .ge. j) then c write(11,*) diffu(j,i) c write(12,*) w(i,j) c else c write(11,*) diffu(i,j) c write(12,*) w(j,i) c endif c enddo c enddo return end