working torsional
[unres4.git] / source / unres / REMD.f90
index 1deb000..8a049a7 100644 (file)
         enddo
         write (iout,*) "Potential forces sidechain"
         do i=nnt,nct
-          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
+          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) &
             write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
         enddo
       endif
        do j=1,3
          ind=1
          do i=nnt,nct
-           if (itype(i).eq.10 .or. itype(i).eq.ntyp1)then
+           if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)then
            rs(ind)=-gcart(j,i)-gxcart(j,i)
             ind=ind+1
           else
 #endif
         ind=1
          do i=nnt,nct
-         if (itype(i).eq.10 .or. itype(i).eq.ntyp1)then 
+         if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)then 
           d_a(j,i)=xsolv(ind)
           ind=ind+1
          else
          d_a(j,0)=d_a(j,nnt)
        enddo
        do i=nnt,nct
-         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
+         if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1) then
            do j=1,3
             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
            enddo
       enddo
       if (lprn) write (iout,*) "Potential forces sidechain"
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
              i,(-gxcart(j,i),j=1,3)
           do j=1,3
         enddo
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
           do j=1,3
             ind=ind+1
             d_a(j,i+nres)=d_a_work(ind)
 !#endif
 !      include 'COMMON.IOUNITS'
 !      include 'COMMON.TIME1'
-      logical :: lprn = .true.
+      logical :: lprn = .false.
       logical :: osob
 #ifndef FIVEDIAG
       real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult
          nind=2
        endif     
        do i=nnt+1,nct-1
-!         if (iabs(itype(i)).eq.ntyp1) cycle
+!         if (iabs(itype(i,1)).eq.ntyp1) cycle
          DM(ind)=2*ip4+mp/2
-         if (iabs(itype(i)).eq.10 .or. iabs(itype(i)).eq.ntyp1) then
-           if (iabs(itype(i)).eq.10) DM(ind)=DM(ind)+msc(10)
+         if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then
+           if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10)
            ind=ind+1
          else
-           DM(ind)=DM(ind)+isc(iabs(itype(i)))
-           DM(ind+1)=msc(iabs(itype(i)))+isc(iabs(itype(i)))
+           DM(ind)=DM(ind)+isc(iabs(itype(i,1)))
+           DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1)))
            ind=ind+2
          endif 
        enddo
        
         ind=1
         do i=nnt,nct
-       if (iabs(itype(i)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
-          DU1(ind)=-isc(iabs(itype(i)))
+       if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
+          DU1(ind)=-isc(iabs(itype(i,1)))
           DU1(ind+1)=0.0d0
          ind=ind+2
        else
 
        ind=1
        do i=nnt,nct-1
-!        if (iabs(itype(i)).eq.ntyp1) cycle
-        write (iout,*) "i",i," itype",itype(i),ntyp1
-       if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
+!        if (iabs(itype(i,1)).eq.ntyp1) cycle
+        write (iout,*) "i",i," itype",itype(i,1),ntyp1
+       if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,1)).ne.ntyp1) then
          DU2(ind)=mp/4-ip4
          DU2(ind+1)=0.0d0
          ind=ind+2
       do i=nnt,nct-1
         ind=ind+1
         ind1=ind1+1
-        coeff=0.25d0*IP
-        massvec(ind1)=mp
+        coeff=0.25d0*IP(1)
+        massvec(ind1)=mp(1)
         Gmat(ind,ind)=coeff
         A(ind1,ind)=0.5d0
       enddo
       m1=nct-nnt+1
       ind=0
       ind1=0
-      msc(ntyp1)=1.0d0
+      msc(ntyp1,1)=1.0d0
       do i=nnt,nct
         ind=ind+1
         ii = ind+m
-        iti=itype(i)
-        massvec(ii)=msc(iabs(iti))
+        iti=itype(i,1)
+        massvec(ii)=msc(iabs(iti),1)
         if (iti.ne.10 .and. iti.ne.ntyp1) then
           ind1=ind1+1
           ii1= ind1+m1
           A(ii,ii1)=1.0d0
-          Gmat(ii1,ii1)=ISC(iabs(iti))
+          Gmat(ii1,ii1)=ISC(iabs(iti),1)
         endif
       enddo
 !  Off-diagonal elements of the dX part of A
       ind=0
       k=nct-nnt
       do i=nnt,nct
-        iti=itype(i)
+        iti=itype(i,1)
         ind=ind+1
         do j=nnt,i
           ii = ind
       Bmat(1,1)=1.0d0
       j=2
       do i=nnt,nct
-        if (itype(i).eq.10) then
+        if (itype(i,1).eq.10) then
           Bmat(i-nnt+2,j-1)=-1
           Bmat(i-nnt+2,j)=1
           j=j+1