Adam's 5D respa update

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

diff --git a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F
index 93fe9ab..0f37efe 100644 (file)
--- a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F
+++ b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F
@@ -81,13 +81,16 @@ C
 c      include 'COMMON.CONTACTS'
       double precision gg(3)
       double precision evdw,evdwij
-      integer i,j,k,itypi,itypj,itypi1,num_conti,iint
+      integer i,j,k,itypi,itypj,itypi1,num_conti,iint,icont
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & sigij,r0ij,rcut,sss1,sssgrad1,sqrij
       double precision sscale,sscagrad
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -97,10 +100,10 @@ c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
+c        do iint=1,nint_gr(i)
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
-          do j=istart(i,iint),iend(i,iint)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -138,8 +141,8 @@ C
                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
               enddo
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
       do i=1,nct
         do j=1,3
@@ -180,13 +183,16 @@ C
 c      include 'COMMON.CONTACTS'
       double precision gg(3)
       double precision evdw,evdwij
-      integer i,j,k,itypi,itypj,itypi1,num_conti,iint
+      integer i,j,k,itypi,itypj,itypi1,num_conti,iint,icont
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
       double precision sscale,sscagrad
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -198,10 +204,10 @@ C Change 12/1/95
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
+c        do iint=1,nint_gr(i)
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
-          do j=istart(i,iint),iend(i,iint)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -235,8 +241,8 @@ C
                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
               enddo
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
       do i=1,nct
         do j=1,3
@@ -274,14 +280,17 @@ C
       include "COMMON.SPLITELE"
       double precision gg(3)
       double precision evdw,evdwij
-      integer i,j,k,itypi,itypj,itypi1,iint
+      integer i,j,k,itypi,itypj,itypi1,iint,icont
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
       double precision sscale,sscagrad
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -291,8 +300,8 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -340,8 +349,8 @@ C
                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
               enddo
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c       enddo        ! iint
       enddo          ! i
       do i=1,nct
         do j=1,3
@@ -370,14 +379,17 @@ C
       include "COMMON.SPLITELE"
       double precision gg(3)
       double precision evdw,evdwij
-      integer i,j,k,itypi,itypj,itypi1,iint
+      integer i,j,k,itypi,itypj,itypi1,iint,icont
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
       double precision sscale,sscagrad
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -387,8 +399,8 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -430,8 +442,8 @@ C
                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
               enddo
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
       do i=1,nct
         do j=1,3
@@ -462,7 +474,7 @@ C
       integer icall
       common /srutu/ icall
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision sss1,sssgrad1
       double precision sscale,sscagrad
@@ -477,7 +489,10 @@ c     else
         lprn=.false.
 c     endif
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -492,8 +507,8 @@ c        dsci_inv=dsc_inv(itypi)
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -559,8 +574,8 @@ C Calculate the angular part of the gradient and sum add the contributions
 C to the appropriate components of the Cartesian gradient.
               call sc_grad_scale((1.0d0-sss)*sss1)
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
 c     stop
       return
@@ -586,7 +601,7 @@ C
       integer icall
       common /srutu/ icall
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision sscale,sscagrad
 c     double precision rrsave(maxdim)
@@ -600,7 +615,10 @@ c     else
         lprn=.false.
 c     endif
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -615,8 +633,8 @@ c        dsci_inv=dsc_inv(itypi)
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -678,8 +696,8 @@ C Calculate the angular part of the gradient and sum add the contributions
 C to the appropriate components of the Cartesian gradient.
               call sc_grad_scale(sss)
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
 c     stop
       return
@@ -706,7 +724,7 @@ C
       logical lprn
       integer xshift,yshift,zshift
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
@@ -720,7 +738,10 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       lprn=.false.
 c     if (icall.eq.0) lprn=.false.
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -743,8 +764,8 @@ c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -898,8 +919,8 @@ C Calculate the radial part of the gradient
 C Calculate angular part of the gradient.
               call sc_grad_scale((1.0d0-sss)*sss1)
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
 c      write (iout,*) "Number of loop steps in EGB:",ind
 cccc      energy_dec=.false.
@@ -927,7 +948,7 @@ C
       logical lprn
       integer xshift,yshift,zshift
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
@@ -941,7 +962,10 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       lprn=.false.
 c     if (icall.eq.0) lprn=.false.
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -964,8 +988,8 @@ c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -1115,8 +1139,8 @@ C Calculate the radial part of the gradient
 C Calculate angular part of the gradient.
               call sc_grad_scale(sss)
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
 c      write (iout,*) "Number of loop steps in EGB:",ind
 cccc      energy_dec=.false.
@@ -1143,7 +1167,7 @@ C
       integer icall
       common /srutu/ icall
       logical lprn
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
      & xi,yi,zi,fac_augm,e_augm
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
@@ -1157,7 +1181,10 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       lprn=.false.
 c     if (icall.eq.0) lprn=.true.
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1172,8 +1199,8 @@ c        dsci_inv=dsc_inv(itypi)
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -1257,8 +1284,8 @@ C Calculate the radial part of the gradient
 C Calculate angular part of the gradient.
               call sc_grad_scale((1.0d0-sss)*sss1)
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
       end
 C-----------------------------------------------------------------------------
@@ -1282,7 +1309,7 @@ C
       integer icall
       common /srutu/ icall
       logical lprn
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
      & xi,yi,zi,fac_augm,e_augm
       double precision evdw
@@ -1296,7 +1323,10 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       lprn=.false.
 c     if (icall.eq.0) lprn=.true.
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1311,8 +1341,8 @@ c        dsci_inv=dsc_inv(itypi)
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -1390,8 +1420,8 @@ C Calculate the radial part of the gradient
 C Calculate angular part of the gradient.
               call sc_grad_scale(sss)
             endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
       end
 C----------------------------------------------------------------------------
@@ -1481,6 +1511,7 @@ C
       include 'COMMON.TIME1'
       include 'COMMON.SHIELD'
       include "COMMON.SPLITELE"
+      integer icont
       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
@@ -1629,7 +1660,10 @@ C     &    .or. itype(i-1).eq.ntyp1
 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
+c      do i=iatel_s,iatel_e
+      do icont=g_listpp_start,g_listpp_end
+        i=newcontlistppi(icont)
+        j=newcontlistppj(icont)
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C     &  .or. itype(i+2).eq.ntyp1
 C     &  .or. itype(i-1).eq.ntyp1
@@ -1653,13 +1687,13 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
 #ifdef FOURBODY
         num_conti=num_cont_hb(i)
 #endif
-        do j=ielstart(i),ielend(i)
+c        do j=ielstart(i),ielend(i)
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
 C     & .or.itype(j+2).eq.ntyp1
 C     & .or.itype(j-1).eq.ntyp1
      &) cycle
           call eelecij_scale(i,j,ees,evdw1,eel_loc)
-        enddo ! j
+c        enddo ! j
 #ifdef FOURBODY
         num_cont_hb(i)=num_conti
 #endif
@@ -2700,12 +2734,16 @@ c      write (iout,*) "evdwpp_short"
       double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
      & dist_temp, dist_init,sss_grad
       double precision sscale,sscagrad
+      integer icont
       evdw1=0.0D0
 C      print *,"WCHODZE"
 c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
 c     & " iatel_e_vdw",iatel_e_vdw
 c      call flush(iout)
-      do i=iatel_s_vdw,iatel_e_vdw
+c      do i=iatel_s_vdw,iatel_e_vdw
+      do icont=g_listpp_vdw_start,g_listpp_vdw_end
+        i=newcontlistpp_vdwi(icont)
+        j=newcontlistpp_vdwj(icont)
         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -2726,7 +2764,7 @@ c      call flush(iout)
 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
 c     &   ' ielend',ielend_vdw(i)
 c        call flush(iout)
-        do j=ielstart_vdw(i),ielend_vdw(i)
+c        do j=ielstart_vdw(i),ielend_vdw(i)
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
           ind=ind+1
           iteli=itel(i)
@@ -2817,7 +2855,7 @@ C            ggg(3)=facvdw*zj
               gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
             enddo
           endif
-        enddo ! j
+c        enddo ! j
       enddo   ! i
       return
       end
@@ -2851,6 +2889,7 @@ C
       double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
      & dist_temp, dist_init
       double precision sscale,sscagrad
+      integer icont
       if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
       evdw2=0.0D0
       evdw2_14=0.0d0
@@ -2859,7 +2898,10 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
 c      if (lprint_short) 
 c     &  write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
 c     & ' iatscp_e=',iatscp_e
-      do i=iatscp_s,iatscp_e
+c      do i=iatscp_s,iatscp_e
+      do icont=g_listscp_start,g_listscp_end
+        i=newcontlistscpi(icont)
+        j=newcontlistscpj(icont)
         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))
@@ -2872,9 +2914,9 @@ c     & ' iatscp_e=',iatscp_e
         zi=mod(zi,boxzsize)
         if (zi.lt.0) zi=zi+boxzsize
 
-        do iint=1,nscp_gr(i)
+c        do iint=1,nscp_gr(i)
 
-        do j=iscpstart(i,iint),iscpend(i,iint)
+c        do j=iscpstart(i,iint),iscpend(i,iint)
           itypj=iabs(itype(j))
           if (itypj.eq.ntyp1) cycle
 C Uncomment following three lines for SC-p interactions
@@ -2969,9 +3011,9 @@ c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
             enddo
           endif
-        enddo
+c        enddo
 
-        enddo ! iint
+c        enddo ! iint
       enddo ! i
       do i=1,nct
         do j=1,3