X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fbrown_step.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fbrown_step.F;h=0be97f546f588114c12dddb23e65a62fbfd3a77d;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git

diff --git a/source/unres/src_MD-M-newcorr/brown_step.F b/source/unres/src_MD-M-newcorr/brown_step.F
new file mode 100644
index 0000000..0be97f5
--- /dev/null
+++ b/source/unres/src_MD-M-newcorr/brown_step.F
@@ -0,0 +1,395 @@
+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
+