Adam's update

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

diff --git a/source/unres/src-HCD-5D/lagrangian_lesyng.F b/source/unres/src-HCD-5D/lagrangian_lesyng.F
index 5de830f..1180645 100644 (file)
--- a/source/unres/src-HCD-5D/lagrangian_lesyng.F
+++ b/source/unres/src-HCD-5D/lagrangian_lesyng.F
@@ -58,7 +58,7 @@ c-------------------------------------------------------------------------
       d_a=0.0d0
       if (lprn) then
         write (iout,*) "Potential forces backbone"
-        do i=nnt,nct
+        do i=1,nres
           write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
         enddo
         write (iout,*) "Potential forces sidechain"
@@ -161,6 +161,8 @@ c      write (iout,*) "Shifting accelerations"
         d_a(:,0)=d_a(:,1)
         d_a(:,1)=0.0d0
       endif
+#define CHUJ
+#ifdef CHUJ
       do ichain=2,nchain
 c        write (iout,*) "ichain",chain_border1(1,ichain)-1,
 c     &     chain_border1(1,ichain)
@@ -175,6 +177,7 @@ c     &   chain_border(2,ichain-1)
      &  d_a(:,chain_border1(1,ichain)-1)+d_a(:,chain_border(2,ichain-1))
         d_a(:,chain_border(2,ichain-1))=0.0d0
       enddo
+#endif
 #else
       inct_prev=0
       do j=1,3
@@ -187,20 +190,6 @@ c     &   chain_border(2,ichain-1)
           d_a(j,inct_prev)=d_a(j,innt)-aaux(j)
         enddo
         inct_prev=inct+1
-#ifdef DEBUG
-        do i=innt,inct
-          if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
-            do j=1,3
-              d_a(j,i)=d_a(j,i+1)-d_a(j,i)
-            enddo
-          else 
-            do j=1,3
-              d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
-              d_a(j,i)=d_a(j,i+1)-d_a(j,i)
-            enddo
-          end if
-        enddo
-#else
         do i=innt,inct
           if (itype(i).ne.10) then
             do j=1,3
@@ -216,17 +205,16 @@ c     &   chain_border(2,ichain-1)
             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
           enddo
         enddo
-#endif
       enddo
 #endif
       if (lprn) then
         write(iout,*) 'acceleration 3D FIVEDIAG in dC and dX'
-        do i=0,nct-1
+        do i=0,nres
           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
         enddo
         do i=nnt,nct
           write (iout,'(i3,3f10.5,3x,3f10.5)') 
-     &     i+nres,(d_a(j,i+nres),j=1,3)
+     &     i,(d_a(j,i+nres),j=1,3)
         enddo
       endif
 #else
@@ -903,10 +891,12 @@ c---------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.LAGRANGE.5diag'
       include 'COMMON.INTERACT'
+      include 'COMMON.VAR'
       integer ndim
       double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
      &  xsolv(ndim),d_a_vec(6*nres)
       integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev
+      accel=0.0d0
       do j=1,3
 Compute accelerations in Calpha and SC
         do ichain=1,nchain
@@ -915,9 +905,18 @@ Compute accelerations in Calpha and SC
           innt=chain_border(1,ichain)
           inct=chain_border(2,ichain)
           do i=iposc,iposc+n-1
-            rs(i)=forces(3*(i-1)+j)
+            rs(i-iposc+1)=forces(3*(i-1)+j)
           enddo
+#ifdef DEBUG
+          write (iout,*) "j",j," chain",ichain
+          write (iout,*) "rs"
+          write (iout,'(f10.5)') (rs(i),i=1,n)
+#endif
           call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
+#ifdef DEBUG
+          write (iout,*) "xsolv"
+          write (iout,'(f10.5)') (xsolv(i),i=1,n)
+#endif
           ind=1
           do i=innt,inct
             if (itype(i).eq.10)then
@@ -931,7 +930,7 @@ Compute accelerations in Calpha and SC
           enddo
         enddo
       enddo
-C Conevert d_a to virtual-bon-vector basis
+C Convert d_a to virtual-bon-vector basis
 #ifdef DEBUG
       write (iout,*) "accel in CA-SC basis"
       do i=1,nres
@@ -956,6 +955,7 @@ C Conevert d_a to virtual-bon-vector basis
         end if
       enddo
       accel(:,nres)=0.0d0
+      accel(:,nct)=0.0d0
       accel(:,2*nres)=0.0d0
       if (nnt.gt.1) then
         accel(:,0)=accel(:,1)
@@ -998,7 +998,7 @@ C Conevert d_a to virtual-bon-vector basis
       enddo
 #ifdef DEBUG
       write (iout,*) "d_a_vec"
-      write (iout,'(3f10.5)') (d_a_vec(j),j=1,dimen3)
+      write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
 #endif
       return
       end