Adam's unres update

[unres.git] / source / unres / src-HCD-5D / stochfric.F

diff --git a/source/unres/src-HCD-5D/stochfric.F b/source/unres/src-HCD-5D/stochfric.F
index dc0b088..946fb58 100644 (file)
--- a/source/unres/src-HCD-5D/stochfric.F
+++ b/source/unres/src-HCD-5D/stochfric.F
@@ -1,5 +1,8 @@
       subroutine friction_force
       implicit none
+#ifdef MPI
+      include 'mpif.h'
+#endif
       include 'DIMENSIONS'
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
@@ -25,7 +28,7 @@
       include 'COMMON.IOUNITS'
 #ifdef FIVEDIAG
       integer iposc,ichain,n,innt,inct
-      double precision v_work(3,maxres2),vvec(maxres2_chain),rs(maxres2)
+      double precision v_work(3,maxres2),vvec(maxres2),rs(maxres2)
 #else
       double precision gamvec(MAXRES6)
       common /syfek/ gamvec
@@ -37,8 +40,12 @@
       
       logical lprn /.false./, checkmode /.false./
 #ifdef FIVEDIAG
+#ifdef TIMING
+      include 'COMMON.TIME1'
+      double precision time01
+#endif
 c Here accelerations due to friction forces are computed right after forces.
-      d_t_work=0.0d0
+      d_t_work(:6*nres)=0.0d0
       do j=1,3
         v_work(j,1)=d_t(j,0)
         v_work(j,nnt)=d_t(j,0)
@@ -70,6 +77,9 @@ c          write (iout,*) "friction_force j",j," ichain",ichain,
 c     &       " n",n," iposc",iposc,iposc+n-1
           innt=chain_border(1,ichain)
           inct=chain_border(2,ichain)
+c diagnostics
+c          innt=chain_border(1,1)
+c          inct=chain_border(2,1)
           do i=innt,inct
             vvec(ind+1)=v_work(j,i)
             ind=ind+1
@@ -79,23 +89,35 @@ c     &       " n",n," iposc",iposc,iposc+n-1
             endif
           enddo
 #ifdef DEBUG
-          write (iout,*) "vvec ind",ind
+          write (iout,*) "vvec ind",ind," n",n
           write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
 #endif
 c          write (iout,*) "chain",i," ind",ind," n",n
+#ifdef TIMING
+          time01=MPI_Wtime()
+#endif
           call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
-     &     DU2fric(iposc),vvec,rs)
+     &     DU2fric(iposc),vvec(iposc),rs)
+#ifdef TIMING
+          time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
+#endif
+#ifdef DEBUG
+          write (iout,*) "rs"
+          write (iout,'(f10.5)') (rs(i),i=1,n)
+#endif
           do i=iposc,iposc+n-1
-            fric_work(3*(i-1)+j)=-rs(i)
+c            write (iout,*) "ichain",ichain," i",i," j",j,
+c     &       "index",3*(i-1)+j,"rs",rs(i-iposc+1)
+            fric_work(3*(i-1)+j)=-rs(i-iposc+1)
           enddo  
         enddo
       enddo
 #ifdef DEBUG
-      write (iout,*) "Vector fric_work"
+      write (iout,*) "Vector fric_work dimen3",dimen3
       write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
 #endif
 #else
-      do i=0,MAXRES2
+      do i=0,2*nres
         do j=1,3
           friction(j,i)=0.0d0
         enddo
@@ -268,7 +290,7 @@ c-----------------------------------------------------
       integer ichain,innt,inct,iposc
 #endif
 
-      do i=0,MAXRES2
+      do i=0,2*nres
         do j=1,3
           stochforc(j,i)=0.0d0
         enddo
@@ -282,6 +304,9 @@ c-----------------------------------------------------
 #endif
 c Compute the stochastic forces acting on bodies. Store in force.
       do i=nnt,nct-1
+#ifdef FIVEDIAG
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+#endif
         sig=stdforcp(i)
         lowb=-5*sig
         highb=5*sig
@@ -315,6 +340,10 @@ c Compute the stochastic forces acting on bodies. Store in force.
         innt=chain_border(1,ichain)
         inct=chain_border(2,ichain)
         iposc=iposd_chain(ichain)
+c for debugging only
+c        innt=chain_border(1,1)
+c        inct=chain_border(2,1)
+c        iposc=iposd_chain(1)
 c        write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
 c     &    " inct",inct," iposc",iposc
         do j=1,3
@@ -417,7 +446,6 @@ c Compute the stochastic forces acting on virtual-bond vectors.
           ind=ind+3
         endif
       enddo
-#endif
       if (lprn) then
         write (iout,*) "stochforcvec"
         do i=1,3*dimen
@@ -468,8 +496,8 @@ c Compute the stochastic forces acting on virtual-bond vectors.
         enddo
         ind=ind+3
       enddo
-
       endif
+#endif
       return
       end
 c------------------------------------------------------------------
@@ -481,6 +509,10 @@ c------------------------------------------------------------------
       double precision time00
 #endif
       include 'DIMENSIONS'
+#ifndef FIVEDIAG
+      integer mmaxres2
+      parameter (mmaxres2=(maxres2*(maxres2+1)/2))
+#endif
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
@@ -542,15 +574,15 @@ C      gamsc(ntyp1)=1.0d0
       enddo
       if (surfarea) call sdarea(gamvec)
       if (lprn) then
-        write (iout,*) "Vector gamvec"
-        do i=1,dimen1
+        write (iout,*) "Vector gamvec ii",ii
+        do i=1,ii
           write (iout,'(i5,f10.5)') i, gamvec(i)
         enddo
       endif
 #ifdef FIVEDIAG
-      DMfric=0.0d0
-      DU1fric=0.0d0
-      DU2fric=0.0d0
+      DMfric(:2*nres)=0.0d0
+      DU1fric(:2*nres)=0.0d0
+      DU2fric(:2*nres)=0.0d0
       ind=1
       do ichain=1,nchain
         innt=chain_border(1,ichain)