zmiany w galezi multichain

[unres.git] / source / unres / src_MD-M / energy_p_new-sep_barrier.F

diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
index 815ca5a..574e19b 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
@@ -13,6 +13,21 @@ C-----------------------------------------------------------------------
       return
       end
 C-----------------------------------------------------------------------
+C-----------------------------------------------------------------------
+      double precision function sscagrad(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+      if(r.lt.r_cut-rlamb) then
+        sscagrad=0.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscagrad=gamm*(6*gamm-6.0d0)/rlamb
+      else
+        sscagrad=0.0d0
+      endif
+      return
+      end
+C-----------------------------------------------------------------------
       subroutine elj_long(evdw)
 C
 C This subroutine calculates the interaction energy of nonbonded side chains
@@ -37,7 +52,7 @@ c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -50,7 +65,7 @@ cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
           do j=istart(i,iint),iend(i,iint)
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -123,7 +138,7 @@ c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -138,7 +153,7 @@ cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
           do j=istart(i,iint),iend(i,iint)
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -209,7 +224,7 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -220,7 +235,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -292,7 +307,7 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -303,7 +318,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -384,7 +399,7 @@ c     endif
       ind=0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -401,7 +416,7 @@ C
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             chi1=chi(itypi,itypj)
@@ -497,7 +512,7 @@ c     endif
       ind=0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -514,7 +529,7 @@ C
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             chi1=chi(itypi,itypj)
@@ -607,7 +622,7 @@ c     if (icall.eq.0) lprn=.false.
       ind=0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -626,7 +641,7 @@ C
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
@@ -745,7 +760,7 @@ c     if (icall.eq.0) lprn=.false.
       ind=0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -764,7 +779,7 @@ C
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
@@ -882,7 +897,7 @@ c     if (icall.eq.0) lprn=.true.
       ind=0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -899,7 +914,7 @@ C
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             sig0ij=sigma(itypi,itypj)
@@ -1004,7 +1019,7 @@ c     if (icall.eq.0) lprn=.true.
       ind=0
       do i=iatsc_s,iatsc_e
         itypi=itype(i)
-        if (itypi.eq.21) cycle
+        if (itypi.eq.ntyp1) cycle
         itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
@@ -1021,7 +1036,7 @@ C
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=itype(j)
-            if (itypj.eq.21) cycle
+            if (itypj.eq.ntyp1) cycle
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             sig0ij=sigma(itypi,itypj)
@@ -1262,8 +1277,8 @@ C
 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
 C
       do i=iturn3_start,iturn3_end
-        if (itype(i).eq.21 .or. itype(i+1).eq.21
-     &  .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle
+        if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
+     &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -1279,9 +1294,9 @@ C
         num_cont_hb(i)=num_conti
       enddo
       do i=iturn4_start,iturn4_end
-        if (itype(i).eq.21 .or. itype(i+1).eq.21
-     &    .or. itype(i+3).eq.21
-     &    .or. itype(i+4).eq.21) cycle
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+     &    .or. itype(i+3).eq.ntyp1
+     &    .or. itype(i+4).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -1293,7 +1308,7 @@ C
         zmedi=c(3,i)+0.5d0*dzi
         num_conti=num_cont_hb(i)
         call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
-        if (wturn4.gt.0.0d0 .and. itype(i+2).ne.21) 
+        if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) 
      &    call eturn4(i,eello_turn4)
         num_cont_hb(i)=num_conti
       enddo   ! i
@@ -1301,7 +1316,7 @@ c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
       do i=iatel_s,iatel_e
-        if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -1314,7 +1329,7 @@ c
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
         do j=ielstart(i),ielend(i)
-          if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle
+          if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
           call eelecij_scale(i,j,ees,evdw1,eel_loc)
         enddo ! j
         num_cont_hb(i)=num_conti
@@ -2000,7 +2015,7 @@ c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
 c     & " iatel_e_vdw",iatel_e_vdw
       call flush(iout)
       do i=iatel_s_vdw,iatel_e_vdw
-        if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
+        if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -2015,7 +2030,7 @@ c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
 c     &   ' ielend',ielend_vdw(i)
         call flush(iout)
         do j=ielstart_vdw(i),ielend_vdw(i)
-          if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle
+          if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
           ind=ind+1
           iteli=itel(i)
           itelj=itel(j)
@@ -2088,7 +2103,7 @@ C
 cd    print '(a)','Enter ESCP'
 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
       do i=iatscp_s,iatscp_e
-        if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
@@ -2098,7 +2113,7 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
 
         do j=iscpstart(i,iint),iscpend(i,iint)
           itypj=itype(j)
-          if (itypj.eq.21) cycle
+          if (itypj.eq.ntyp1) cycle
 C Uncomment following three lines for SC-p interactions
 c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
@@ -2189,7 +2204,7 @@ C
 cd    print '(a)','Enter ESCP'
 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
       do i=iatscp_s,iatscp_e
-        if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
@@ -2199,7 +2214,7 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
 
         do j=iscpstart(i,iint),iscpend(i,iint)
           itypj=itype(j)
-          if (itypj.eq.21) cycle
+          if (itypj.eq.ntyp1) cycle
 C Uncomment following three lines for SC-p interactions
 c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi