X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=77103ccdfd340a0cbd75f8b25ba0151a9af52bbf;hb=9978233a8f362a38c9d5929a421ac612aed6bf7e;hp=7960906c5e540e9a88a4c8222297d8bd1e249a9d;hpb=b05d82c3d028d22db2d3b9bb47e40f0c59be753f;p=unres4.git

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index 7960906..77103cc 100644
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -180,7 +180,7 @@
       subroutine etotal(energia)
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
-      use MD_data, only: totT
+      use MD_data
 #ifndef ISNAN
       external proc_proc
 #ifdef WINPGI
@@ -213,6 +213,14 @@
 
 #ifdef MPI      
       real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
+! shielding effect varibles for MPI
+!      real(kind=8)   fac_shieldbuf(maxres),
+!     & grad_shield_locbuf(3,maxcontsshi,-1:maxres),
+!     & grad_shield_sidebuf(3,maxcontsshi,-1:maxres),
+!     & grad_shieldbuf(3,-1:maxres)
+!       integer ishield_listbuf(maxres),
+!     &shield_listbuf(maxcontsshi,maxres)
+
 !      print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
 !     & " nfgtasks",nfgtasks
       if (nfgtasks.gt.1) then
@@ -1269,9 +1277,11 @@
 !      include 'COMMON.SBRIDGE'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj
+      integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
       real(kind=8) :: evdw,sig0ij
+      real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+                    dist_temp, dist_init
       integer :: ii
 !cccc      energy_dec=.false.
 !     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -1286,6 +1296,13 @@
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+          xi=dmod(xi,boxxsize)
+          if (xi.lt.0) xi=xi+boxxsize
+          yi=dmod(yi,boxysize)
+          if (yi.lt.0) yi=yi+boxysize
+          zi=dmod(zi,boxzsize)
+          if (zi.lt.0) zi=zi+boxzsize
+
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1334,9 +1351,46 @@
 !           alf1=0.0D0
 !           alf2=0.0D0
 !           alf12=0.0D0
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+           xj=c(1,nres+j)
+           yj=c(2,nres+j)
+           zj=c(3,nres+j)
+          xj=dmod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=dmod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=dmod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      subchap=0
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            subchap=1
+          endif
+       enddo
+       enddo
+       enddo
+       if (subchap.eq.1) then
+          xj=xj_temp-xi
+          yj=yj_temp-yi
+          zj=zj_temp-zi
+       else
+          xj=xj_safe-xi
+          yj=yj_safe-yi
+          zj=zj_safe-zi
+       endif
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -1349,6 +1403,11 @@
 !          write(iout,*)"c ", c(1,:), c(2,:), c(3,:)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
+            sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+!            print *,sss_ele_cut,sss_ele_grad,&
+!            1.0d0/(rij),r_cut_ele,rlamb_ele
+            if (sss_ele_cut.le.0.0) cycle
 ! Calculate angle-dependent terms of energy and contributions to their
 ! derivatives.
             call sc_angular
@@ -1380,7 +1439,7 @@
 !          write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,& !d
 !          " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2," fac",fac !d
             evdwij=evdwij*eps2rt*eps3rt
-            evdw=evdw+evdwij
+            evdw=evdw+evdwij*sss_ele_cut
             if (lprn) then
             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -1402,11 +1461,18 @@
             fac=-expon*(e1+evdwij)*rij_shift
             sigder=fac*sigder
             fac=rij*fac
+!            print *,'before fac',fac,rij,evdwij
+            fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij
+!            print *,'grad part scale',fac,   &
+!             evdwij*sss_ele_grad/sss_ele_cut &
+!            /sigma(itypi,itypj)*rij
 !            fac=0.0d0
 ! Calculate the radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
+!            print *,'before sc_grad', gg(1),gg(2),gg(3)
 ! Calculate angular part of the gradient.
             call sc_grad
             ENDIF    ! dyn_ss            
@@ -2625,6 +2691,12 @@
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=dmod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=dmod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=dmod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
         num_conti=0
         call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
@@ -2643,6 +2715,12 @@
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=dmod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=dmod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=dmod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
         num_conti=num_cont_hb(i)
         call eelecij(i,i+3,ees,evdw1,eel_loc)
         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) &
@@ -2663,6 +2741,13 @@
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=dmod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=dmod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=dmod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+
 !        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
         do j=ielstart(i),ielend(i)
@@ -2709,6 +2794,9 @@
       real(kind=8),dimension(2,2) :: acipa !el,a_temp
 !el      real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
       real(kind=8),dimension(4) :: muij
+      real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+                    dist_temp, dist_init
+      integer xshift,yshift,zshift
 !el      integer :: num_conti,j1,j2
 !el      real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
 !el        dz_normi,xmedi,ymedi,zmedi
@@ -2730,7 +2818,7 @@
                                              0.0d0,0.0d0,1.0d0/),shape(unmat)) 
 !      integer :: maxconts=nres/4
 !el local variables
-      integer :: k,i,j,iteli,itelj,kkk,l,kkll,m
+      integer :: k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap
       real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
       real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i
       real(kind=8) :: dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,&
@@ -2761,12 +2849,54 @@
           dx_normj=dc_norm(1,j)
           dy_normj=dc_norm(2,j)
           dz_normj=dc_norm(3,j)
-          xj=c(1,j)+0.5D0*dxj-xmedi
-          yj=c(2,j)+0.5D0*dyj-ymedi
-          zj=c(3,j)+0.5D0*dzj-zmedi
+!          xj=c(1,j)+0.5D0*dxj-xmedi
+!          yj=c(2,j)+0.5D0*dyj-ymedi
+!          zj=c(3,j)+0.5D0*dzj-zmedi
+          xj=c(1,j)+0.5D0*dxj
+          yj=c(2,j)+0.5D0*dyj
+          zj=c(3,j)+0.5D0*dzj
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            isubchap=1
+          endif
+       enddo
+       enddo
+       enddo
+       if (isubchap.eq.1) then
+!C          print *,i,j
+          xj=xj_temp-xmedi
+          yj=yj_temp-ymedi
+          zj=zj_temp-zmedi
+       else
+          xj=xj_safe-xmedi
+          yj=yj_safe-ymedi
+          zj=zj_safe-zmedi
+       endif
+
           rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
+            sss_ele_cut=sscale_ele(rij)
+            sss_ele_grad=sscagrad_ele(rij)
+!            print *,sss_ele_cut,sss_ele_grad,&
+!            1.0d0/(rij),r_cut_ele,rlamb_ele
+            if (sss_ele_cut.le.0.0) go to 128
+
           rmij=1.0D0/rij
           r3ij=rrmij*rmij
           r6ij=r3ij*r3ij  
@@ -2786,8 +2916,8 @@
           eesij=el1+el2
 ! 12/26/95 - for the evaluation of multi-body H-bonding interactions
           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
-          ees=ees+eesij
-          evdw1=evdw1+evdwij
+          ees=ees+eesij*sss_ele_cut
+          evdw1=evdw1+evdwij*sss_ele_cut
 !d          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 !d     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
 !d     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
@@ -2804,8 +2934,8 @@
 ! Calculate contributions to the Cartesian gradient.
 !
 #ifdef SPLITELE
-          facvdw=-6*rrmij*(ev1+evdwij)
-          facel=-3*rrmij*(el1+eesij)
+          facvdw=-6*rrmij*(ev1+evdwij)*sss_ele_cut
+          facel=-3*rrmij*(el1+eesij)*sss_ele_cut
           fac1=fac
           erij(1)=xj*rmij
           erij(2)=yj*rmij
@@ -2834,9 +2964,9 @@
 !grad              gelc(l,k)=gelc(l,k)+ggg(l)
 !grad            enddo
 !grad          enddo
-          ggg(1)=facvdw*xj
-          ggg(2)=facvdw*yj
-          ggg(3)=facvdw*zj
+          ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj
+          ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj
+          ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj
 !          do k=1,3
 !            ghalf=0.5D0*ggg(k)
 !            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
@@ -2856,8 +2986,8 @@
 !grad            enddo
 !grad          enddo
 #else
-          facvdw=ev1+evdwij 
-          facel=el1+eesij  
+          facvdw=(ev1+evdwij)*sss_ele_cut
+          facel=(el1+eesij)*sss_ele_cut
           fac1=fac
           fac=-3*rrmij*(facvdw+facvdw+facel)
           erij(1)=xj*rmij
@@ -2866,9 +2996,9 @@
 !
 ! Radial derivatives. First process both termini of the fragment (i,j)
 ! 
-          ggg(1)=fac*xj
-          ggg(2)=fac*yj
-          ggg(3)=fac*zj
+          ggg(1)=fac*xj+sss_ele_grad*rmij*(eesij+evdwij)*xj
+          ggg(2)=fac*yj+sss_ele_grad*rmij*(eesij+evdwij)*yj
+          ggg(3)=fac*zj+sss_ele_grad*rmij*(eesij+evdwij)*zj
 !          do k=1,3
 !            ghalf=0.5D0*ggg(k)
 !            gelc(k,i)=gelc(k,i)+ghalf
@@ -2937,6 +3067,7 @@
             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
           enddo
+ 128      continue
           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 &
               .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 &
               .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
@@ -4119,16 +4250,19 @@
 !      if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres)
 
       do i=ibondp_start,ibondp_end
+        if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
         if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
-          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
-          do j=1,3
-          gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
-            *dc(j,i-1)/vbld(i)
-          enddo
-          if (energy_dec) write(iout,*) &
-             "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+!C          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+!C          do j=1,3
+!C          gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
+!C            *dc(j,i-1)/vbld(i)
+!C          enddo
+!C          if (energy_dec) write(iout,*) &
+!C             "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+        diff = vbld(i)-vbldpDUM
         else
         diff = vbld(i)-vbldp0
+        endif
         if (energy_dec) write (iout,'(a7,i5,4f7.3)') &
            "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
         estr=estr+diff*diff
@@ -4136,7 +4270,7 @@
           gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
         enddo
 !        write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
-        endif
+!        endif
       enddo
       estr=0.5d0*AKP*estr+estr1
 !
@@ -9564,21 +9698,27 @@
 !       " sigder",sigder
 !      write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
 !      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
+!C      print *,sss_ele_cut,'in sc_grad'
       do k=1,3
         dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
       enddo
       do k=1,3
-        gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+        gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss_ele_cut
+!C      print *,'gg',k,gg(k)
       enddo 
 !      write (iout,*) "gg",(gg(k),k=1,3)
       do k=1,3
         gvdwx(k,i)=gvdwx(k,i)-gg(k) &
                   +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
-                  +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+                  +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv    &
+                  *sss_ele_cut
+
         gvdwx(k,j)=gvdwx(k,j)+gg(k) &
                   +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
-                  +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+                  +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv    &
+                  *sss_ele_cut
+
 !        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
 !                 +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 !        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
@@ -10228,6 +10368,257 @@
       enddo
       return
       end subroutine check_ecart
+#ifdef CARGRAD
+!-----------------------------------------------------------------------------
+      subroutine check_ecartint
+! Check the gradient of the energy in Cartesian coordinates. 
+      use io_base, only: intout
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.CONTROL'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.DERIV'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.VAR'
+!      include 'COMMON.CONTACTS'
+!      include 'COMMON.MD'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.SPLITELE'
+      use comm_srutu
+!el      integer :: icall
+!el      common /srutu/ icall
+      real(kind=8),dimension(6) :: ggg,ggg1
+      real(kind=8),dimension(3) :: cc,xx,ddc,ddx,ddc1,ddcn
+      real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
+      real(kind=8),dimension(3) :: dcnorm_safe1,dcnorm_safe2,dxnorm_safe
+      real(kind=8),dimension(6,0:nres) :: grad_s,grad_s1 !(6,0:maxres)
+      real(kind=8),dimension(nres) :: phi_temp,theta_temp,alph_temp,omeg_temp !(maxres)
+      real(kind=8),dimension(0:n_ene) :: energia,energia1
+      integer :: uiparm(1)
+      real(kind=8) :: urparm(1)
+!EL      external fdum
+      integer :: i,j,k,nf
+      real(kind=8) :: rlambd,aincr,etot,etot1,etot11,etot12,etot2,&
+                   etot21,etot22
+      r_cut=2.0d0
+      rlambd=0.3d0
+      icg=1
+      nf=0
+      nfl=0
+      call intout
+!      call intcartderiv
+!      call checkintcartgrad
+      call zerograd
+      aincr=1.0D-5
+      write(iout,*) 'Calling CHECK_ECARTINT.'
+      nf=0
+      icall=0
+      write (iout,*) "Before geom_to_var"
+      call geom_to_var(nvar,x)
+      write (iout,*) "after geom_to_var"
+      write (iout,*) "split_ene ",split_ene
+      call flush(iout)
+      if (.not.split_ene) then
+        write(iout,*) 'Calling CHECK_ECARTINT if'
+        call etotal(energia)
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+        etot=energia(0)
+        write (iout,*) "etot",etot
+        call flush(iout)
+!el        call enerprint(energia)
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+        call flush(iout)
+        write (iout,*) "enter cartgrad"
+        call flush(iout)
+        call cartgrad
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+        write (iout,*) "exit cartgrad"
+        call flush(iout)
+        icall =1
+        do i=1,nres
+          write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
+        enddo
+        do j=1,3
+          grad_s(j,0)=gcart(j,0)
+        enddo
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+        do i=1,nres
+          do j=1,3
+            grad_s(j,i)=gcart(j,i)
+            grad_s(j+3,i)=gxcart(j,i)
+          enddo
+        enddo
+      else
+write(iout,*) 'Calling CHECK_ECARTIN else.'
+!- split gradient check
+        call zerograd
+        call etotal_long(energia)
+!el        call enerprint(energia)
+        call flush(iout)
+        write (iout,*) "enter cartgrad"
+        call flush(iout)
+        call cartgrad
+        write (iout,*) "exit cartgrad"
+        call flush(iout)
+        icall =1
+        write (iout,*) "longrange grad"
+        do i=1,nres
+          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+          (gxcart(j,i),j=1,3)
+        enddo
+        do j=1,3
+          grad_s(j,0)=gcart(j,0)
+        enddo
+        do i=1,nres
+          do j=1,3
+            grad_s(j,i)=gcart(j,i)
+            grad_s(j+3,i)=gxcart(j,i)
+          enddo
+        enddo
+        call zerograd
+        call etotal_short(energia)
+!el        call enerprint(energia)
+        call flush(iout)
+        write (iout,*) "enter cartgrad"
+        call flush(iout)
+        call cartgrad
+        write (iout,*) "exit cartgrad"
+        call flush(iout)
+        icall =1
+        write (iout,*) "shortrange grad"
+        do i=1,nres
+          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+          (gxcart(j,i),j=1,3)
+        enddo
+        do j=1,3
+          grad_s1(j,0)=gcart(j,0)
+        enddo
+        do i=1,nres
+          do j=1,3
+            grad_s1(j,i)=gcart(j,i)
+            grad_s1(j+3,i)=gxcart(j,i)
+          enddo
+        enddo
+      endif
+      write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
+!      do i=1,nres
+      do i=nnt,nct
+        do j=1,3
+          if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1)
+          if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres)
+	  ddc(j)=c(j,i) 
+	  ddx(j)=c(j,i+nres) 
+          dcnorm_safe1(j)=dc_norm(j,i-1)
+          dcnorm_safe2(j)=dc_norm(j,i)
+          dxnorm_safe(j)=dc_norm(j,i+nres)
+        enddo
+	do j=1,3
+	  c(j,i)=ddc(j)+aincr
+          if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr
+          if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr
+          if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          call int_from_cart1(.false.)
+          if (.not.split_ene) then
+            call etotal(energia1)
+            etot1=energia1(0)
+            write (iout,*) "ij",i,j," etot1",etot1
+          else
+!- split gradient
+            call etotal_long(energia1)
+            etot11=energia1(0)
+            call etotal_short(energia1)
+            etot12=energia1(0)
+          endif
+!- end split gradient
+!          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
+	  c(j,i)=ddc(j)-aincr
+          if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr
+          if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr
+          if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          call int_from_cart1(.false.)
+          if (.not.split_ene) then
+            call etotal(energia1)
+            etot2=energia1(0)
+            write (iout,*) "ij",i,j," etot2",etot2
+	    ggg(j)=(etot1-etot2)/(2*aincr)
+          else
+!- split gradient
+            call etotal_long(energia1)
+            etot21=energia1(0)
+	    ggg(j)=(etot11-etot21)/(2*aincr)
+            call etotal_short(energia1)
+            etot22=energia1(0)
+	    ggg1(j)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+!            write (iout,*) "etot21",etot21," etot22",etot22
+          endif
+!          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+	  c(j,i)=ddc(j)
+          if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)
+          if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)
+          if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          dc_norm(j,i-1)=dcnorm_safe1(j)
+          dc_norm(j,i)=dcnorm_safe2(j)
+          dc_norm(j,i+nres)=dxnorm_safe(j)
+        enddo
+	do j=1,3
+	  c(j,i+nres)=ddx(j)+aincr
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          call int_from_cart1(.false.)
+          if (.not.split_ene) then
+            call etotal(energia1)
+            etot1=energia1(0)
+          else
+!- split gradient
+            call etotal_long(energia1)
+            etot11=energia1(0)
+            call etotal_short(energia1)
+            etot12=energia1(0)
+          endif
+!- end split gradient
+	  c(j,i+nres)=ddx(j)-aincr
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          call int_from_cart1(.false.)
+          if (.not.split_ene) then
+            call etotal(energia1)
+            etot2=energia1(0)
+	    ggg(j+3)=(etot1-etot2)/(2*aincr)
+          else
+!- split gradient
+            call etotal_long(energia1)
+            etot21=energia1(0)
+	    ggg(j+3)=(etot11-etot21)/(2*aincr)
+            call etotal_short(energia1)
+            etot22=energia1(0)
+	    ggg1(j+3)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+          endif
+!          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+	  c(j,i+nres)=ddx(j)
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          dc_norm(j,i+nres)=dxnorm_safe(j)
+          call int_from_cart1(.false.)
+        enddo
+	write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+         i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
+        if (split_ene) then
+          write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+         i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),&
+         k=1,6)
+         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+         i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),&
+         ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
+        endif
+      enddo
+      return
+      end subroutine check_ecartint
+#else
 !-----------------------------------------------------------------------------
       subroutine check_ecartint
 ! Check the gradient of the energy in Cartesian coordinates. 
@@ -10268,7 +10659,7 @@
 !      call intcartderiv
 !      call checkintcartgrad
       call zerograd
-      aincr=1.0D-5
+      aincr=1.0D-6
       write(iout,*) 'Calling CHECK_ECARTINT.'
       nf=0
       icall=0
@@ -10466,6 +10857,7 @@
       enddo
       return
       end subroutine check_ecartint
+#endif
 !-----------------------------------------------------------------------------
       subroutine check_eint
 ! Check the gradient of energy in internal coordinates.
@@ -10672,6 +11064,49 @@
       endif
       return
       end function sscale
+      real(kind=8) function sscale_grad(r)
+!      include "COMMON.SPLITELE"
+      real(kind=8) :: r,gamm
+      if(r.lt.r_cut-rlamb) then
+        sscale_grad=0.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscale_grad=gamm*(6*gamm-6.0d0)/rlamb
+      else
+        sscale_grad=0d0
+      endif
+      return
+      end function sscale_grad
+
+!!!!!!!!!! PBCSCALE
+      real(kind=8) function sscale_ele(r)
+!      include "COMMON.SPLITELE"
+      real(kind=8) :: r,gamm
+      if(r.lt.r_cut_ele-rlamb_ele) then
+        sscale_ele=1.0d0
+      else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+        gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+        sscale_ele=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+      else
+        sscale_ele=0d0
+      endif
+      return
+      end function sscale_ele
+
+      real(kind=8)  function sscagrad_ele(r)
+      real(kind=8) :: r,gamm
+!      include "COMMON.SPLITELE"
+      if(r.lt.r_cut_ele-rlamb_ele) then
+        sscagrad_ele=0.0d0
+      else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+        gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+        sscagrad_ele=gamm*(6*gamm-6.0d0)/rlamb_ele
+      else
+        sscagrad_ele=0.0d0
+      endif
+      return
+      end function sscagrad_ele
+!!!!!!!!!!!!!!!
 !-----------------------------------------------------------------------------
       subroutine elj_long(evdw)
 !
@@ -11290,9 +11725,12 @@
 !      include 'COMMON.CONTROL'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj
+      integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift
-      real(kind=8) :: sss,e1,e2,evdw
+      real(kind=8) :: sss,e1,e2,evdw,sss_grad
+      real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+                    dist_temp, dist_init
+
       evdw=0.0D0
 !cccc      energy_dec=.false.
 !     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -11307,6 +11745,12 @@
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+          xi=mod(xi,boxxsize)
+          if (xi.lt.0) xi=xi+boxxsize
+          yi=mod(yi,boxysize)
+          if (yi.lt.0) yi=yi+boxysize
+          zi=mod(zi,boxzsize)
+          if (zi.lt.0) zi=zi+boxzsize
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -11337,16 +11781,58 @@
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+! Searching for nearest neighbour
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+          dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          xj_safe=xj
+          yj_safe=yj
+          zj_safe=zj
+          subchap=0
+          do xshift=-1,1
+          do yshift=-1,1
+          do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            subchap=1
+          endif
+          enddo
+          enddo
+          enddo
+          if (subchap.eq.1) then
+          xj=xj_temp-xi
+          yj=yj_temp-yi
+          zj=zj_temp-zi
+          else
+          xj=xj_safe-xi
+          yj=yj_safe-yi
+          zj=zj_safe-zi
+          endif
+
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
-
+            sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
+            if (sss_ele_cut.le.0.0) cycle
             if (sss.lt.1.0d0) then
 
 ! Calculate angle-dependent terms of energy and contributions to their
@@ -11377,7 +11863,7 @@
 !              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 !     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
               evdwij=evdwij*eps2rt*eps3rt
-              evdw=evdw+evdwij*(1.0d0-sss)
+              evdw=evdw+evdwij*(1.0d0-sss)*sss_ele_cut
               if (lprn) then
               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -11399,6 +11885,9 @@
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
+              fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij  &
+            /sigmaii(itypi,itypj))
 !              fac=0.0d0
 ! Calculate the radial part of the gradient
               gg(1)=xj*fac
@@ -11435,9 +11924,11 @@
 !      include 'COMMON.CONTROL'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj
+      integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
-      real(kind=8) :: sss,e1,e2,evdw,rij_shift
+      real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad
+      real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+                    dist_temp, dist_init
       evdw=0.0D0
 !cccc      energy_dec=.false.
 !     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -11452,6 +11943,18 @@
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+          xi=mod(xi,boxxsize)
+          if (xi.lt.0) xi=xi+boxxsize
+          yi=mod(yi,boxysize)
+          if (yi.lt.0) yi=yi+boxysize
+          zi=mod(zi,boxzsize)
+          if (zi.lt.0) zi=zi+boxzsize
+        dxi=dc_norm(1,nres+i)
+        dyi=dc_norm(2,nres+i)
+        dzi=dc_norm(3,nres+i)
+!        dsci_inv=dsc_inv(itypi)
+        dsci_inv=vbld_inv(i+nres)
+
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -11482,15 +11985,61 @@
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+!            xj=c(1,nres+j)-xi
+!            yj=c(2,nres+j)-yi
+!            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+! Searching for nearest neighbour
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+          dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          xj_safe=xj
+          yj_safe=yj
+          zj_safe=zj
+          subchap=0
+          do xshift=-1,1
+          do yshift=-1,1
+          do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            subchap=1
+          endif
+          enddo
+          enddo
+          enddo
+          if (subchap.eq.1) then
+          xj=xj_temp-xi
+          yj=yj_temp-yi
+          zj=zj_temp-zi
+          else
+          xj=xj_safe-xi
+          yj=yj_safe-yi
+          zj=zj_safe-zi
+          endif
+
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+            sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
+            sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            if (sss_ele_cut.le.0.0) cycle
 
             if (sss.gt.0.0d0) then
 
@@ -11522,7 +12071,7 @@
 !              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 !     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
               evdwij=evdwij*eps2rt*eps3rt
-              evdw=evdw+evdwij*sss
+              evdw=evdw+evdwij*sss*sss_ele_cut
               if (lprn) then
               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -11544,6 +12093,10 @@
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
+              fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij+sss_grad/sss*rij  &
+            /sigmaii(itypi,itypj))
+
 !              fac=0.0d0
 ! Calculate the radial part of the gradient
               gg(1)=xj*fac
@@ -13039,16 +13592,19 @@
         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
       enddo
       do k=1,3
-        gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
+        gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac&
+         *sss_ele_cut
       enddo 
 !      write (iout,*) "gg",(gg(k),k=1,3)
       do k=1,3
         gvdwx(k,i)=gvdwx(k,i)-gg(k) &
                   +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
-                +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
+                +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac&
+                 *sss_ele_cut
         gvdwx(k,j)=gvdwx(k,j)+gg(k) &
                   +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
-                +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
+                +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac&
+         *sss_ele_cut
 !        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
 !     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 !        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
@@ -13072,7 +13628,7 @@
 !
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
-      use MD_data, only: totT
+      use MD_data, only: totT,usampl,eq_time
 #ifndef ISNAN
       external proc_proc
 #ifdef WINPGI
@@ -13752,7 +14308,7 @@
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
       use energy_data
-      use MD_data, only: totT
+      use MD_data, only: totT,usampl,eq_time
 #ifdef MPI
       include 'mpif.h'
 #endif
@@ -13839,6 +14395,32 @@
             (gxcart(j,i),j=1,3)
       enddo
 #endif
+#ifdef CARGRAD
+#ifdef DEBUG
+      write (iout,*) "CARGRAD"
+#endif
+      do i=nres,1,-1
+        do j=1,3
+          gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+!          gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+        enddo
+!        write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
+!            (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
+      enddo    
+! Correction: dummy residues
+        if (nnt.gt.1) then
+          do j=1,3
+!            gcart_new(j,nnt)=gcart_new(j,nnt)+gcart_new(j,1)
+            gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
+          enddo
+        endif
+        if (nct.lt.nres) then
+          do j=1,3
+!            gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
+            gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
+          enddo
+        endif
+#endif
 #ifdef TIMING
       time_cartgrad=time_cartgrad+MPI_Wtime()-time00
 #endif
@@ -15896,7 +16478,7 @@
 !-----------------------------------------------------------------------------
       subroutine alloc_ener_arrays
 !EL Allocation of arrays used by module energy
-
+      use MD_data, only: mset
 !el local variables
       integer :: i,j
       
@@ -15983,10 +16565,8 @@
       allocate(mu(2,nres))
       allocate(muder(2,nres))
       allocate(Ub2(2,nres))
-        do i=1,nres
-          Ub2(1,i)=0.0d0
-          Ub2(2,i)=0.0d0
-        enddo
+      Ub2(1,:)=0.0d0
+      Ub2(2,:)=0.0d0
       allocate(Ub2der(2,nres))
       allocate(Ctobr(2,nres))
       allocate(Ctobrder(2,nres))
@@ -16169,9 +16749,7 @@
 !      allocate(qinfrag(50,nprocs/20),wfrag(50,nprocs/20)) !(50,maxprocs/20)
 !      allocate(qinpair(100,nprocs/20),wpair(100,nprocs/20)) !(100,maxprocs/20)
       allocate(mset(0:nprocs))  !(maxprocs/20)
-      do i=0,nprocs
-        mset(i)=0
-      enddo
+      mset(:)=0
 !      allocate(ifrag(2,50,nprocs/20))  !(2,50,maxprocs/20)
 !      allocate(ipair(2,100,nprocs/20))  !(2,100,maxprocs/20)
       allocate(dUdconst(3,0:nres))
@@ -16190,11 +16768,11 @@
 ! and side-chain vectors in theta or phi.
       allocate(dyn_ssbond_ij(0:nres+4,0:nres+4))
 !(maxres,maxres)
-      do i=1,nres
-        do j=i+1,nres
-          dyn_ssbond_ij(i,j)=1.0d300
-        enddo
-      enddo
+!      do i=1,nres
+!        do j=i+1,nres
+      dyn_ssbond_ij(:,:)=1.0d300
+!        enddo
+!      enddo
 
       if (nss.gt.0) then
         allocate(idssb(nss),jdssb(nss))
@@ -16202,9 +16780,7 @@
       endif
       allocate(dyn_ss_mask(nres))
 !(maxres)
-      do i=1,nres
-        dyn_ss_mask(i)=.false.
-      enddo
+      dyn_ss_mask(:)=.false.
 !----------------------
 ! common.sccor
 ! Parameters of the SCCOR term