--- /dev/null
+c-------------------------------------------------------------------------------
+ subroutine brown_step(itime)
+c------------------------------------------------
+c Perform a single Euler integration step of Brownian dynamics
+c------------------------------------------------
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ double precision zapas(MAXRES6)
+ integer ilen,rstcount
+ external ilen
+ double precision stochforcvec(MAXRES6)
+ double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
+ & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
+ & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
+ & ppvec(maxres2)
+ common /stochcalc/ stochforcvec
+ common /gucio/ cm
+ integer itime
+ logical lprn /.false./,lprn1 /.false./
+ integer maxiter /5/
+ double precision difftol /1.0d-5/
+ nbond=nct-nnt
+ do i=nnt,nct
+ if (itype(i).ne.10) nbond=nbond+1
+ enddo
+c
+ if (lprn1) then
+ write (iout,*) "Generalized inverse of fricmat"
+ call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
+ endif
+ do i=1,dimen
+ do j=1,nbond
+ Bmat(i,j)=0.0d0
+ enddo
+ enddo
+ ind=3
+ ind1=0
+ do i=nnt,nct-1
+ ind1=ind1+1
+ do j=1,3
+ Bmat(ind+j,ind1)=dC_norm(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind1=ind1+1
+ do j=1,3
+ Bmat(ind+j,ind1)=dC_norm(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Bmat"
+ call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
+ endif
+ do i=1,dimen
+ do j=1,nbond
+ GBmat(i,j)=0.0d0
+ do k=1,dimen
+ GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix GBmat"
+ call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
+ endif
+ do i=1,nbond
+ do j=1,nbond
+ Cmat(i,j)=0.0d0
+ do k=1,dimen
+ Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Cmat"
+ call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
+ endif
+ call matinvert(nbond,MAXRES2,Cmat,Cinv)
+ if (lprn1) then
+ write (iout,*) "Matrix Cinv"
+ call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
+ endif
+ do i=1,dimen
+ do j=1,nbond
+ Tmat(i,j)=0.0d0
+ do k=1,nbond
+ Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Tmat"
+ call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
+ endif
+ do i=1,dimen
+ do j=1,dimen
+ if (i.eq.j) then
+ Pmat(i,j)=1.0d0
+ else
+ Pmat(i,j)=0.0d0
+ endif
+ do k=1,nbond
+ Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Pmat"
+ call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
+ endif
+ do i=1,dimen
+ Td(i)=0.0d0
+ ind=0
+ do k=nnt,nct-1
+ ind=ind+1
+ Td(i)=Td(i)+vbl*Tmat(i,ind)
+ enddo
+ do k=nnt,nct
+ if (itype(k).ne.10) then
+ ind=ind+1
+ Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
+ endif
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Vector Td"
+ do i=1,dimen
+ write (iout,'(i5,f10.5)') i,Td(i)
+ enddo
+ endif
+ call stochastic_force(stochforcvec)
+ if (lprn) then
+ write (iout,*) "stochforcvec"
+ do i=1,dimen
+ write (iout,*) i,stochforcvec(i)
+ enddo
+ endif
+ do j=1,3
+ zapas(j)=-gcart(j,0)+stochforcvec(j)
+ d_t_work(j)=d_t(j,0)
+ dC_work(j)=dC_old(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ ind=ind+1
+ zapas(ind)=-gcart(j,i)+stochforcvec(ind)
+ dC_work(ind)=dC_old(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ ind=ind+1
+ zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
+ dC_work(ind)=dC_old(j,i+nres)
+ enddo
+ endif
+ enddo
+
+ if (lprn) then
+ write (iout,*) "Initial d_t_work"
+ do i=1,dimen
+ write (iout,*) i,d_t_work(i)
+ enddo
+ endif
+
+ do i=1,dimen
+ d_t_work(i)=0.0d0
+ do j=1,dimen
+ d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
+ enddo
+ enddo
+
+ do i=1,dimen
+ zapas(i)=Td(i)
+ do j=1,dimen
+ zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Final d_t_work and zapas"
+ do i=1,dimen
+ write (iout,*) i,d_t_work(i),zapas(i)
+ enddo
+ endif
+
+ do j=1,3
+ d_t(j,0)=d_t_work(j)
+ dc(j,0)=zapas(j)
+ dc_work(j)=dc(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_work(i)
+ dc(j,i)=zapas(ind+j)
+ dc_work(ind+j)=dc(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ do j=1,3
+ d_t(j,i+nres)=d_t_work(ind+j)
+ dc(j,i+nres)=zapas(ind+j)
+ dc_work(ind+j)=dc(j,i+nres)
+ enddo
+ ind=ind+3
+ enddo
+ if (lprn) then
+ call chainbuild_cart
+ write (iout,*) "Before correction for rotational lengthening"
+ write (iout,*) "New coordinates",
+ & " and differences between actual and standard bond lengths"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ xx=vbld(i+1)-vbl
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
+ & i,(dC(j,i),j=1,3),xx
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ xx=vbld(i+nres)-vbldsc0(1,itype(i))
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
+ & i,(dC(j,i+nres),j=1,3),xx
+ endif
+ enddo
+ endif
+c Second correction (rotational lengthening)
+c do iter=1,maxiter
+ diffmax=0.0d0
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ blen2 = scalar(dc(1,i),dc(1,i))
+ ppvec(ind)=2*vbl**2-blen2
+ diffbond=dabs(vbl-dsqrt(blen2))
+ if (diffbond.gt.diffmax) diffmax=diffbond
+ if (ppvec(ind).gt.0.0d0) then
+ ppvec(ind)=dsqrt(ppvec(ind))
+ else
+ ppvec(ind)=0.0d0
+ endif
+ if (lprn) then
+ write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
+ endif
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
+ ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
+ diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
+ if (diffbond.gt.diffmax) diffmax=diffbond
+ if (ppvec(ind).gt.0.0d0) then
+ ppvec(ind)=dsqrt(ppvec(ind))
+ else
+ ppvec(ind)=0.0d0
+ endif
+ if (lprn) then
+ write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
+ endif
+ endif
+ enddo
+ if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
+ if (diffmax.lt.difftol) goto 10
+ do i=1,dimen
+ Td(i)=0.0d0
+ do j=1,nbond
+ Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
+ enddo
+ enddo
+ do i=1,dimen
+ zapas(i)=Td(i)
+ do j=1,dimen
+ zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
+ enddo
+ enddo
+ do j=1,3
+ dc(j,0)=zapas(j)
+ dc_work(j)=zapas(j)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ dc(j,i)=zapas(ind+j)
+ dc_work(ind+j)=zapas(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ dc(j,i+nres)=zapas(ind+j)
+ dc_work(ind+j)=zapas(ind+j)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+c Building the chain from the newly calculated coordinates
+ call chainbuild_cart
+ if(ntwe.ne.0) then
+ if (large.and. mod(itime,ntwe).eq.0) then
+ write (iout,*) "Cartesian and internal coordinates: step 1"
+ call cartprint
+ call intout
+ write (iout,'(a)') "Potential forces"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
+ & (-gxcart(j,i),j=1,3)
+ enddo
+ write (iout,'(a)') "Stochastic forces"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
+ & (stochforc(j,i+nres),j=1,3)
+ enddo
+ write (iout,'(a)') "Velocities"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
+ & (d_t(j,i+nres),j=1,3)
+ enddo
+ endif
+ endif
+ if (lprn) then
+ write (iout,*) "After correction for rotational lengthening"
+ write (iout,*) "New coordinates",
+ & " and differences between actual and standard bond lengths"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ xx=vbld(i+1)-vbl
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
+ & i,(dC(j,i),j=1,3),xx
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ xx=vbld(i+nres)-vbldsc0(1,itype(i))
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
+ & i,(dC(j,i+nres),j=1,3),xx
+ endif
+ enddo
+ endif
+c ENDDO
+c write (iout,*) "Too many attempts at correcting the bonds"
+c stop
+ 10 continue
+#ifdef MPI
+ tt0 =MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+c Calculate energy and forces
+ call zerograd
+ call etotal(potEcomp)
+ potE=potEcomp(0)-potEcomp(20)
+ call cartgrad
+ totT=totT+d_time
+c Calculate the kinetic and total energy and the kinetic temperature
+ call kinetic(EK)
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+ totE=EK+potE
+ kinetic_T=2.0d0/(dimen*Rb)*EK
+ return
+ end
+
projects / unres.git / blobdiff