X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=a7b07987ebba738b947d71905a22d3cdacf914a2;hp=44e5cba1dbc6e633cce1abf56814388f34e960f4;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hpb=3a74889c3f6514588d40ec326a78f33b3663f5be

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index 44e5cba..a7b0798 100644
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -44,11 +44,13 @@ C Gay-Berne potential (shifted LJ, angular dependence).
       goto 106
 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
   105 call egbv(evdw,evdw_t)
+C      write(iout,*) 'po elektostatyce'
 C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
-  106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-C
+  106  call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+C            write(iout,*) 'po eelec'
+
 C Calculate excluded-volume interaction energy between peptide groups
 C and side chains.
 C
@@ -56,8 +58,9 @@ C
 c
 c Calculate the bond-stretching energy
 c
+
       call ebond(estr)
-c      write (iout,*) "estr",estr
+C       write (iout,*) "estr",estr
 C 
 C Calculate the disulfide-bridge and other energy and the contributions
 C from other distance constraints.
@@ -67,13 +70,14 @@ cd    print *,'EHPB exitted succesfully.'
 C
 C Calculate the virtual-bond-angle energy.
 C
+C      print *,'Bend energy finished.'
       call ebend(ebe,ethetacnstr)
 cd    print *,'Bend energy finished.'
 C
 C Calculate the SC local energy.
 C
       call esc(escloc)
-cd    print *,'SCLOC energy finished.'
+C       print *,'SCLOC energy finished.'
 C
 C Calculate the virtual-bond torsional energy.
 C
@@ -811,6 +815,14 @@ c      if (icall.gt.0) lprn=.true.
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+C returning the ith atom to box
+          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)
@@ -864,15 +876,59 @@ c           chip12=0.0D0
 c           alf1=0.0D0
 c           alf2=0.0D0
 c           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)
+C returning jth atom to box
+          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
+C checking the distance
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      subchap=0
+C finding the closest
+      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)
 c            write (iout,*) i,j,xj,yj,zj
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
+            sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
+            sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
+            if (sss.le.0.0) cycle
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
             call sc_angular
@@ -895,9 +951,9 @@ c---------------------------------------------------------------
             eps3der=evdwij*eps2rt
             evdwij=evdwij*eps2rt*eps3rt
             if (bb(itypi,itypj).gt.0) then
-              evdw=evdw+evdwij
+              evdw=evdw+evdwij*sss
             else
-              evdw_t=evdw_t+evdwij
+              evdw_t=evdw_t+evdwij*sss
             endif
             ij=icant(itypi,itypj)
             aux=eps1*eps2rt**2*eps3rt**2
@@ -926,6 +982,7 @@ C Calculate gradient components.
             fac=-expon*(e1+evdwij)*rij_shift
             sigder=fac*sigder
             fac=rij*fac
+            fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
 C Calculate the radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
@@ -1881,14 +1938,22 @@ cd      enddo
       do i=1,nres
         num_cont_hb(i)=0
       enddo
-cd      print '(a)','Enter EELEC'
-cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
+C      print '(a)','Enter EELEC'
+C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
       do i=1,nres
         gel_loc_loc(i)=0.0d0
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+          if (i.eq.1) then 
+           if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
+     &  .or. itype(i+2).eq.ntyp1) cycle
+          else
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+     &  .or. itype(i+2).eq.ntyp1
+     &  .or. itype(i-1).eq.ntyp1
+     &) cycle
+         endif
         if (itel(i).eq.0) goto 1215
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -1899,10 +1964,27 @@ cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=mod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=mod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=mod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
         num_conti=0
-c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+C        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
-          if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+          if (j.eq.1) then
+           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+     & .or.itype(j+2).eq.ntyp1
+     &) cycle  
+          else     
+          if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+     & .or.itype(j+2).eq.ntyp1
+     & .or.itype(j-1).eq.ntyp1
+     &) cycle
+         endif
+C
+C) cycle
           if (itel(j).eq.0) goto 1216
           ind=ind+1
           iteli=itel(i)
@@ -1924,10 +2006,49 @@ C End diagnostics
           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
+          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)
+          if (zj.lt.0) zj=zj+boxzsize
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      isubchap=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-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
+          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
+            sss=sscale(sqrt(rij))
+            sssgrad=sscagrad(sqrt(rij))
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
           rmij=1.0D0/rij
@@ -1951,12 +2072,12 @@ c          write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij
 C 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
+          evdw1=evdw1+evdwij*sss
 c             write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') 
 c     &'evdw1',i,j,evdwij
 c     &,iteli,itelj,aaa,evdw1
 
-c              write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+C              write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
 c          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 c     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
 c     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
@@ -1965,7 +2086,7 @@ C
 C Calculate contributions to the Cartesian gradient.
 C
 #ifdef SPLITELE
-          facvdw=-6*rrmij*(ev1+evdwij) 
+          facvdw=-6*rrmij*(ev1+evdwij)*sss
           facel=-3*rrmij*(el1+eesij)
           fac1=fac
           erij(1)=xj*rmij
@@ -1991,9 +2112,18 @@ C
               gelc(l,k)=gelc(l,k)+ggg(l)
             enddo
           enddo
-          ggg(1)=facvdw*xj
-          ggg(2)=facvdw*yj
-          ggg(3)=facvdw*zj
+C          ggg(1)=facvdw*xj
+C          ggg(2)=facvdw*yj
+C          ggg(3)=facvdw*zj
+          if (sss.gt.0.0) then
+          ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+          ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+          ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+          else
+          ggg(1)=0.0
+          ggg(2)=0.0
+          ggg(3)=0.0
+          endif
           do k=1,3
             ghalf=0.5D0*ggg(k)
             gvdwpp(k,i)=gvdwpp(k,i)+ghalf
@@ -2008,7 +2138,7 @@ C
             enddo
           enddo
 #else
-          facvdw=ev1+evdwij 
+          facvdw=(ev1+evdwij)*sss
           facel=el1+eesij  
           fac1=fac
           fac=-3*rrmij*(facvdw+facvdw+facel)
@@ -2854,7 +2984,13 @@ c     &   " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-
+C Returning the ith atom to box
+          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
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
@@ -2865,28 +3001,73 @@ c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
 c         zj=c(3,nres+j)-zi
 C Uncomment following three lines for Ca-p interactions
-          xj=c(1,j)-xi
-          yj=c(2,j)-yi
-          zj=c(3,j)-zi
+          xj=c(1,j)
+          yj=c(2,j)
+          zj=c(3,j)
+C returning the jth atom to box
+          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
+C Finding the closest jth atom
+      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
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+C sss is scaling function for smoothing the cutoff gradient otherwise
+C the gradient would not be continuouse
+          sss=sscale(1.0d0/(dsqrt(rrij)))
+          if (sss.le.0.0d0) cycle
+          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
           fac=rrij**expon2
           e1=fac*fac*aad(itypj,iteli)
           e2=fac*bad(itypj,iteli)
           if (iabs(j-i) .le. 2) then
             e1=scal14*e1
             e2=scal14*e2
-            evdw2_14=evdw2_14+e1+e2
+            evdw2_14=evdw2_14+(e1+e2)*sss
           endif
           evdwij=e1+e2
 c          write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
 c     &        'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
 c     &       bad(itypj,iteli)
-          evdw2=evdw2+evdwij
+          evdw2=evdw2+evdwij*sss
           if (calc_grad) then
 C
 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
 C
-          fac=-(evdwij+e1)*rrij
+          fac=-(evdwij+e1)*rrij*sss
+          fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
           ggg(1)=xj*fac
           ggg(2)=yj*fac
           ggg(3)=zj*fac
@@ -3197,24 +3378,29 @@ c
       estr1=0.0d0
 c      write (iout,*) "distchainmax",distchainmax
       do i=nnt+1,nct
-        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,vbld(i),distchainmax,
-     &       gnmr1(vbld(i),-1.0d0,distchainmax)
-        else
+        if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
+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,vbld(i),distchainmax,
+C     &       gnmr1(vbld(i),-1.0d0,distchainmax)
+C        else
+         if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
+        diff = vbld(i)-vbldpDUM
+         else
           diff = vbld(i)-vbldp0
 c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
+         endif
           estr=estr+diff*diff
           do j=1,3
             gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
           enddo
-        endif
-
+C        endif
+C        write (iout,'(a7,i5,4f7.3)')
+C     &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
       enddo
       estr=0.5d0*AKP*estr+estr1
 c
@@ -3226,8 +3412,8 @@ c
           nbi=nbondterm(iti)
           if (nbi.eq.1) then
             diff=vbld(i+nres)-vbldsc0(1,iti)
-c            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
-c     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
+C            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+C     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
             estr=estr+0.5d0*AKSC(1,iti)*diff*diff
             do j=1,3
               gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
@@ -3299,7 +3485,10 @@ c      write (iout,*) "nres",nres
 c     write (*,'(a,i2)') 'EBEND ICG=',icg
 c      write (iout,*) ithet_start,ithet_end
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
+C        if (itype(i-1).eq.ntyp1) cycle
+        if (i.le.2) cycle
+        if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+     &  .or.itype(i).eq.ntyp1) cycle
 C Zero the energy function and its derivative at 0 or pi.
         call splinthet(theta(i),0.5d0*delta,ss,ssd)
         it=itype(i-1)
@@ -3315,8 +3504,12 @@ C Zero the energy function and its derivative at 0 or pi.
           ichir21=isign(1,itype(i))
           ichir22=isign(1,itype(i))
          endif
+         if (i.eq.3) then
+          y(1)=0.0D0
+          y(2)=0.0D0
+          else
 
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
 c          icrc=0
@@ -3331,7 +3524,8 @@ c          call proc_proc(phii,icrc)
           y(1)=0.0D0
           y(2)=0.0D0
         endif
-        if (i.lt.nres .and. itype(i).ne.ntyp1) then
+        endif
+        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
 c          icrc=0
@@ -3398,6 +3592,8 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
      &        E_theta,E_tc)
         endif
         etheta=etheta+ethetai
+c         write (iout,'(a6,i5,0pf7.3,f7.3,i5)')
+c     &      'ebend',i,ethetai,theta(i),itype(i)
 c        write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i),
 c     &    rad2deg*phii,rad2deg*phii1,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
@@ -3574,9 +3770,11 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
-c        if (itype(i-1).eq.ntyp1) cycle
-        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
-     &(itype(i).eq.ntyp1)) cycle
+C         if (i.eq.2) cycle
+C        if (itype(i-1).eq.ntyp1) cycle
+        if (i.le.2) cycle
+        if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+     &  .or.itype(i).eq.ntyp1) cycle
         if (iabs(itype(i+1)).eq.20) iblock=2
         if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
@@ -3588,6 +3786,14 @@ c        if (itype(i-1).eq.ntyp1) cycle
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
+        if (i.eq.3) then 
+          phii=0.0d0
+          ityp1=nthetyp+1
+          do k=1,nsingle
+            cosph1(k)=0.0d0
+            sinph1(k)=0.0d0
+          enddo
+        else
         if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
@@ -3609,6 +3815,7 @@ c          ityp1=nthetyp+1
             sinph1(k)=0.0d0
           enddo 
         endif
+        endif
         if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
@@ -3799,14 +4006,14 @@ C ALPHA and OMEGA.
       common /sccalc/ time11,time12,time112,theti,it,nlobit
       delta=0.02d0*pi
       escloc=0.0D0
-c     write (iout,'(a)') 'ESC'
+C      write (iout,*) 'ESC'
       do i=loc_start,loc_end
         it=itype(i)
         if (it.eq.ntyp1) cycle
         if (it.eq.10) goto 1
         nlobit=nlob(iabs(it))
 c       print *,'i=',i,' it=',it,' nlobit=',nlobit
-c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
+C        write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
         theti=theta(i+1)-pipol
         x(1)=dtan(theti)
         x(2)=alph(i)
@@ -3842,8 +4049,8 @@ c        write (iout,*) "i",i," x",x(1),x(2),x(3)
             dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
           enddo
           dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
-c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c    &             esclocbi,ss,ssd
+          write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
+     &             esclocbi,ss,ssd
           escloci=ss*escloci+(1.0d0-ss)*esclocbi
 c         escloci=esclocbi
 c         write (iout,*) escloci
@@ -3877,15 +4084,17 @@ c         write (iout,*) escloci
           enddo
           dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
 c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c    &             esclocbi,ss,ssd
+c     &             esclocbi,ss,ssd
           escloci=ss*escloci+(1.0d0-ss)*esclocbi
-c         write (iout,*) escloci
+C         write (iout,*) 'i=',i, escloci
         else
           call enesc(x,escloci,dersc,ddummy,.false.)
         endif
 
         escloc=escloc+escloci
-c        write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
+C        write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
+            write (iout,'(a6,i5,0pf7.3)')
+     &     'escloc',i,escloci
 
         gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
      &   wscloc*dersc(1)
@@ -4576,8 +4785,11 @@ C Set lprn=.true. for debugging
 c      lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
-        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
-     &       .or. itype(i).eq.ntyp1) cycle
+        if (i.le.2) cycle
+        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+     &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+C        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
+C     &       .or. itype(i).eq.ntyp1) cycle
         if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215
          if (iabs(itype(i)).eq.20) then
          iblock=2
@@ -4680,8 +4892,12 @@ C Set lprn=.true. for debugging
 c     lprn=.true.
       etors_d=0.0D0
       do i=iphi_start,iphi_end-1
-        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
-     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (i.le.3) cycle
+C        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+C     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+         if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or.
+     &  (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or.
+     &  (itype(i+1).eq.ntyp1)) cycle
         if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0) 
      &     goto 1215
         itori=itortyp(itype(i-2))
@@ -7691,4 +7907,34 @@ C-----------------------------------------------------------------------------
       scalar=sc
       return
       end
+C-----------------------------------------------------------------------
+      double precision function sscale(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+      if(r.lt.r_cut-rlamb) then
+        sscale=1.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+      else
+        sscale=0d0
+      endif
+      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-----------------------------------------------------------------------