adding ebend_nucl to UCGM+some further reading

[unres4.git] / source / unres / energy.f90
diff --git a/source/unres/energy.f90 b/source/unres/energy.f90

index c166edf..2fc2773 100644 (file)
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -45,6 +45,7 @@
        real(kind=8),dimension(:,:,:),allocatable :: gacont      !(3,maxconts,maxres)
        integer,dimension(:),allocatable :: ishield_list
        integer,dimension(:,:),allocatable ::  shield_list
        real(kind=8),dimension(:,:,:),allocatable :: gacont      !(3,maxconts,maxres)
        integer,dimension(:),allocatable :: ishield_list
        integer,dimension(:,:),allocatable ::  shield_list
+      real(kind=8),dimension(:),allocatable :: enetube,enecavtube
  !                
  ! 12/26/95 - H-bonding contacts
  !      common /contacts_hb/ 
  !                
  ! 12/26/95 - H-bonding contacts
  !      common /contacts_hb/ 
@@ -123,7 +124,9 @@
          gshieldc_ec,gshieldc_loc_ec,gshieldx_t3, &
          gshieldc_t3,gshieldc_loc_t3,gshieldx_t4,gshieldc_t4, &
          gshieldc_loc_t4,gshieldx_ll,gshieldc_ll,gshieldc_loc_ll,&
          gshieldc_ec,gshieldc_loc_ec,gshieldx_t3, &
          gshieldc_t3,gshieldc_loc_t3,gshieldx_t4,gshieldc_t4, &
          gshieldc_loc_t4,gshieldx_ll,gshieldc_ll,gshieldc_loc_ll,&
-        grad_shield,gg_tube,gg_tube_sc !(3,maxres)
+        grad_shield,gg_tube,gg_tube_sc,gradafm !(3,maxres)
+!-----------------------------NUCLEIC GRADIENT
+      real(kind=8),dimension(:,:),allocatable  ::gradb_nucl,gradbx_nucl
  !      real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
        real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,&
          gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres)
  !      real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
        real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,&
          gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres)
@@ -222,9 +225,13 @@
        integer :: n_corr,n_corr1,ierror
        real(kind=8) :: etors,edihcnstr,etors_d,esccor,ehpb
        real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc
        integer :: n_corr,n_corr1,ierror
        real(kind=8) :: etors,edihcnstr,etors_d,esccor,ehpb
        real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc
-      real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube
+      real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, &
+                      Eafmforce,ethetacnstr
        real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6
        real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6
-
+! now energies for nulceic alone parameters
+      real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
+                      ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
+                      ecorr3_nucl
  #ifdef MPI      
        real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
  ! shielding effect varibles for MPI
  #ifdef MPI      
        real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
  ! shielding effect varibles for MPI
@@ -307,7 +314,7 @@
  #endif
  ! 
  ! Compute the side-chain and electrostatic interaction energy
  #endif
  ! 
  ! Compute the side-chain and electrostatic interaction energy
-        print *, "Before EVDW"
+!        print *, "Before EVDW"
  !      goto (101,102,103,104,105,106) ipot
        select case(ipot)
  ! Lennard-Jones potential.
  !      goto (101,102,103,104,105,106) ipot
        select case(ipot)
  ! Lennard-Jones potential.
@@ -355,6 +362,7 @@
         if (shield_mode.eq.2) then
                   call set_shield_fac2
         endif
         if (shield_mode.eq.2) then
                   call set_shield_fac2
         endif
+       print *,"AFTER EGB",ipot,evdw
  !mc
  !mc Sep-06: egb takes care of dynamic ss bonds too
  !mc
  !mc
  !mc Sep-06: egb takes care of dynamic ss bonds too
  !mc
@@ -370,7 +378,7 @@
  !        print *,"Processor",myrank," left VEC_AND_DERIV"
        if (ipot.lt.6) then
  #ifdef SPLITELE
  !        print *,"Processor",myrank," left VEC_AND_DERIV"
        if (ipot.lt.6) then
  #ifdef SPLITELE
-         print *,"after ipot if", ipot
+!         print *,"after ipot if", ipot
           if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or. &
               wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0 &
               .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
           if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or. &
               wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0 &
               .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
@@ -420,6 +428,7 @@
  ! Calculate the bond-stretching energy
  !
        call ebond(estr)
  ! Calculate the bond-stretching energy
  !
        call ebond(estr)
+       print *,"EBOND",estr
  !       write(iout,*) "in etotal afer ebond",ipot
  
  ! 
  !       write(iout,*) "in etotal afer ebond",ipot
  
  ! 
@@ -433,7 +442,7 @@
  ! Calculate the virtual-bond-angle energy.
  !
        if (wang.gt.0d0) then
  ! Calculate the virtual-bond-angle energy.
  !
        if (wang.gt.0d0) then
-        call ebend(ebe)
+        call ebend(ebe,ethetacnstr)
        else
          ebe=0
        endif
        else
          ebe=0
        endif
@@ -520,6 +529,13 @@
        else
         eliptran=0.0d0
        endif
        else
         eliptran=0.0d0
        endif
+      if (fg_rank.eq.0) then
+      if (AFMlog.gt.0) then
+        call AFMforce(Eafmforce)
+      else if (selfguide.gt.0) then
+        call AFMvel(Eafmforce)
+      endif
+      endif
        if (tubemode.eq.1) then
         call calctube(etube)
        else if (tubemode.eq.2) then
        if (tubemode.eq.1) then
         call calctube(etube)
        else if (tubemode.eq.2) then
@@ -529,7 +545,10 @@
        else
         etube=0.0d0
        endif
        else
         etube=0.0d0
        endif
-
+!--------------------------------------------------------
+      call ebond_nucl(estr_nucl)
+      call ebend_nucl(ebe_nucl)
+      print *,"after ebend", ebe_nucl
  #ifdef TIMING
        time_enecalc=time_enecalc+MPI_Wtime()-time00
  #endif
  #ifdef TIMING
        time_enecalc=time_enecalc+MPI_Wtime()-time00
  #endif
@@ -572,7 +591,24 @@
        energia(20)=Uconst+Uconst_back
        energia(21)=esccor
        energia(22)=eliptran
        energia(20)=Uconst+Uconst_back
        energia(21)=esccor
        energia(22)=eliptran
+      energia(23)=Eafmforce
+      energia(24)=ethetacnstr
        energia(25)=etube
        energia(25)=etube
+!---------------------------------------------------------------
+      energia(26)=evdwpp
+      energia(27)=eespp
+      energia(28)=evdwpsb
+      energia(29)=eelpsb
+      energia(30)=evdwsb
+      energia(31)=eelsb
+      energia(32)=estr_nucl
+      energia(33)=ebe_nucl
+      energia(34)=esbloc
+      energia(35)=etors_nucl
+      energia(36)=etors_d_nucl
+      energia(37)=ecorr_nucl
+      energia(38)=ecorr3_nucl
+!----------------------------------------------------------------------
  !    Here are the energies showed per procesor if the are more processors 
  !    per molecule then we sum it up in sum_energy subroutine 
  !      print *," Processor",myrank," calls SUM_ENERGY"
  !    Here are the energies showed per procesor if the are more processors 
  !    per molecule then we sum it up in sum_energy subroutine 
  !      print *," Processor",myrank," calls SUM_ENERGY"
@@ -614,7 +650,11 @@
        real(kind=8) :: evdw,evdw2,evdw2_14,ees,evdw1,ecorr,ecorr5,ecorr6
        real(kind=8) :: eel_loc,eello_turn3,eello_turn4,eturn6,ebe,escloc
        real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot,   &
        real(kind=8) :: evdw,evdw2,evdw2_14,ees,evdw1,ecorr,ecorr5,ecorr6
        real(kind=8) :: eel_loc,eello_turn3,eello_turn4,eturn6,ebe,escloc
        real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot,   &
-        eliptran,etube
+        eliptran,etube, Eafmforce,ethetacnstr
+      real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
+                      ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
+                      ecorr3_nucl
+
        integer :: i
  #ifdef MPI
        integer :: ierr
        integer :: i
  #ifdef MPI
        integer :: ierr
@@ -675,21 +715,30 @@
        Uconst=energia(20)
        esccor=energia(21)
        eliptran=energia(22)
        Uconst=energia(20)
        esccor=energia(21)
        eliptran=energia(22)
+      Eafmforce=energia(23)
+      ethetacnstr=energia(24)
        etube=energia(25)
        etube=energia(25)
+      estr_nucl=energia(32)
+      ebe_nucl=energia(33)
+
  #ifdef SPLITELE
        etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 &
         +wang*ebe+wtor*etors+wscloc*escloc &
         +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
         +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 &
         +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
  #ifdef SPLITELE
        etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 &
         +wang*ebe+wtor*etors+wscloc*escloc &
         +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
         +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 &
         +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
-       +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube
+       +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+       +Eafmforce+ethetacnstr  &
+       +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl
  #else
        etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) &
         +wang*ebe+wtor*etors+wscloc*escloc &
         +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
         +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 &
         +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
  #else
        etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) &
         +wang*ebe+wtor*etors+wscloc*escloc &
         +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
         +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 &
         +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
-       +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube
+       +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+       +Eafmforce+ethetacnstr &
+       +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl
  #endif
        energia(0)=etot
  ! detecting NaNQ
  #endif
        energia(0)=etot
  ! detecting NaNQ
@@ -813,7 +862,10 @@
        real(kind=8) :: etot,evdw,evdw2,ees,evdw1,ecorr,ecorr5,ecorr6,eel_loc
        real(kind=8) :: eello_turn6,eello_turn3,eello_turn4,ebe,escloc
        real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran,&
        real(kind=8) :: etot,evdw,evdw2,ees,evdw1,ecorr,ecorr5,ecorr6,eel_loc
        real(kind=8) :: eello_turn6,eello_turn3,eello_turn4,ebe,escloc
        real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran,&
-       etube
+       etube,ethetacnstr,Eafmforce
+      real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
+                      ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
+                      ecorr3_nucl
  
        etot=energia(0)
        evdw=energia(1)
  
        etot=energia(0)
        evdw=energia(1)
@@ -844,7 +896,12 @@
        Uconst=energia(20)
        esccor=energia(21)
        eliptran=energia(22)
        Uconst=energia(20)
        esccor=energia(21)
        eliptran=energia(22)
+      Eafmforce=energia(23)
+      ethetacnstr=energia(24)
        etube=energia(25)
        etube=energia(25)
+      estr_nucl=energia(32)
+      ebe_nucl=energia(33)
+
  #ifdef SPLITELE
        write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
          estr,wbond,ebe,wang,&
  #ifdef SPLITELE
        write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
          estr,wbond,ebe,wang,&
@@ -852,8 +909,10 @@
          ecorr,wcorr,&
          ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,&
          eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,&
          ecorr,wcorr,&
          ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,&
          eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,&
-        edihcnstr,ebr*nss,&
-        Uconst,eliptran,wliptran,etube,wtube,etot
+        edihcnstr,ethetacnstr,ebr*nss,&
+        Uconst,eliptran,wliptran,Eafmforce,etube,wtube, & ! till now protein
+        estr_nucl,wbond_nucl,ebe_nucl,wang_nucl, &
+        etot
     10 format (/'Virtual-chain energies:'// &
         'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
         'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
     10 format (/'Virtual-chain energies:'// &
         'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
         'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
@@ -875,10 +934,14 @@
         'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ &
         'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ &
         'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ &
         'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ &
         'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ &
         'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ &
+       'ETHETC= ',1pE16.6,' (valence angle constraints)'/ &
         'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ &
         'UCONST= ',1pE16.6,' (Constraint energy)'/ &
         'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/&
         'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ &
         'UCONST= ',1pE16.6,' (Constraint energy)'/ &
         'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/&
+       'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/ &
         'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ &
         'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ &
+       'ESTR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (stretching for nucleic)'/ &
+       'EBE_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (bending for nucleic)'/ &
         'ETOT=  ',1pE16.6,' (total)')
  #else
        write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,&
         'ETOT=  ',1pE16.6,' (total)')
  #else
        write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,&
@@ -887,7 +950,10 @@
          ecorr,wcorr,&
          ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,&
          eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,&
          ecorr,wcorr,&
          ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,&
          eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,&
-        ebr*nss,Uconst,eliptran,wliptran,etube,wtube,etot
+        ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,     &
+        etube,wtube, &
+        estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,&
+        etot
     10 format (/'Virtual-chain energies:'// &
         'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
         'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
     10 format (/'Virtual-chain energies:'// &
         'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
         'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
@@ -908,10 +974,14 @@
         'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ &
         'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ &
         'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ &
         'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ &
         'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ &
         'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ &
+       'ETHETC= ',1pE16.6,' (valence angle constraints)'/ &
         'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ &
         'UCONST=',1pE16.6,' (Constraint energy)'/ &
         'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ &
         'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ &
         'UCONST=',1pE16.6,' (Constraint energy)'/ &
         'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ &
+       'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/ &
         'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ &
         'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ &
+       'ESTR_nucl=  ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching for nucleic)'/ &
+       'EBE_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (bending for nucleic)'/ &
         'ETOT=  ',1pE16.6,' (total)')
  #endif
        return
         'ETOT=  ',1pE16.6,' (total)')
  #endif
        return
@@ -952,9 +1022,9 @@
  !      allocate(gacont(3,nres/4,iatsc_s:iatsc_e))      !(3,maxconts,maxres)
  
        do i=iatsc_s,iatsc_e
  !      allocate(gacont(3,nres/4,iatsc_s:iatsc_e))      !(3,maxconts,maxres)
  
        do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
+        itypi=iabs(itype(i,1))
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=iabs(itype(i+1))
+        itypi1=iabs(itype(i+1,1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -967,7 +1037,7 @@
  !d        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  !d   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
  !d        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  !d   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
-            itypj=iabs(itype(j)) 
+            itypj=iabs(itype(j,1)) 
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
@@ -984,7 +1054,7 @@
  !d          sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
  !d          epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
  !d          write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
  !d          sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
  !d          epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
  !d          write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d   &        restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d   &        restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
  !d   &        bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
  !d   &        (c(k,i),k=1,3),(c(k,j),k=1,3)
              evdw=evdw+evdwij
  !d   &        bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
  !d   &        (c(k,i),k=1,3),(c(k,j),k=1,3)
              evdw=evdw+evdwij
@@ -1110,9 +1180,9 @@
  !     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
  !     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
+        itypi=iabs(itype(i,1))
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=iabs(itype(i+1))
+        itypi1=iabs(itype(i+1,1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -1121,7 +1191,7 @@
  !
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
-            itypj=iabs(itype(j))
+            itypj=iabs(itype(j,1))
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
@@ -1139,7 +1209,7 @@
  !d          sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
  !d          epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
  !d          write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
  !d          sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
  !d          epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
  !d          write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d   &        restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d   &        restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
  !d   &        bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
  !d   &        sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
  !d   &        (c(k,i),k=1,3),(c(k,j),k=1,3)
  !d   &        bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
  !d   &        sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
  !d   &        (c(k,i),k=1,3),(c(k,j),k=1,3)
@@ -1211,9 +1281,9 @@
  !     endif
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !     endif
  !el      ind=0
        do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
+        itypi=iabs(itype(i,1))
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=iabs(itype(i+1))
+        itypi1=iabs(itype(i+1,1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -1228,7 +1298,7 @@
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
-            itypj=iabs(itype(j))
+            itypj=iabs(itype(j,1))
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
@@ -1280,7 +1350,7 @@
              sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
              epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
  !d            write (iout,'(2(a3,i3,2x),15(0pf7.3))')
              sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
              epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
  !d            write (iout,'(2(a3,i3,2x),15(0pf7.3))')
-!d     &        restyp(itypi),i,restyp(itypj),j,
+!d     &        restyp(itypi,1),i,restyp(itypj,1),j,
  !d     &        epsi,sigm,chi1,chi2,chip1,chip2,
  !d     &        eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
  !d     &        om1,om2,om12,1.0D0/dsqrt(rrij),
  !d     &        epsi,sigm,chi1,chi2,chip1,chip2,
  !d     &        eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
  !d     &        om1,om2,om12,1.0D0/dsqrt(rrij),
@@ -1343,10 +1413,10 @@
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !C        print *,"I am in EVDW",i
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !C        print *,"I am in EVDW",i
-        itypi=iabs(itype(i))
+        itypi=iabs(itype(i,1))
  !        if (i.ne.47) cycle
          if (itypi.eq.ntyp1) cycle
  !        if (i.ne.47) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=iabs(itype(i+1))
+        itypi1=iabs(itype(i+1,1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -1379,7 +1449,7 @@
           sslipi=0.0d0
           ssgradlipi=0.0
         endif
           sslipi=0.0d0
           ssgradlipi=0.0
         endif
-       print *, sslipi,ssgradlipi
+!       print *, sslipi,ssgradlipi
          dxi=dc_norm(1,nres+i)
          dyi=dc_norm(2,nres+i)
          dzi=dc_norm(3,nres+i)
          dxi=dc_norm(1,nres+i)
          dyi=dc_norm(2,nres+i)
          dzi=dc_norm(3,nres+i)
@@ -1399,16 +1469,37 @@
                                'evdw',i,j,evdwij,' ss'
  !              if (energy_dec) write (iout,*) &
  !                              'evdw',i,j,evdwij,' ss'
                                'evdw',i,j,evdwij,' ss'
  !              if (energy_dec) write (iout,*) &
  !                              'evdw',i,j,evdwij,' ss'
+             do k=j+1,iend(i,iint)
+!C search over all next residues
+              if (dyn_ss_mask(k)) then
+!C check if they are cysteins
+!C              write(iout,*) 'k=',k
+
+!c              write(iout,*) "PRZED TRI", evdwij
+!               evdwij_przed_tri=evdwij
+              call triple_ssbond_ene(i,j,k,evdwij)
+!c               if(evdwij_przed_tri.ne.evdwij) then
+!c                 write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+!c               endif
+
+!c              write(iout,*) "PO TRI", evdwij
+!C call the energy function that removes the artifical triple disulfide
+!C bond the soubroutine is located in ssMD.F
+              evdw=evdw+evdwij
+              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+                            'evdw',i,j,evdwij,'tss'
+              endif!dyn_ss_mask(k)
+             enddo! k
              ELSE
  !el            ind=ind+1
              ELSE
  !el            ind=ind+1
-            itypj=iabs(itype(j))
+            itypj=iabs(itype(j,1))
              if (itypj.eq.ntyp1) cycle
  !             if (j.ne.78) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
  !            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,&
  !              1.0d0/vbld(j+nres) !d
              if (itypj.eq.ntyp1) cycle
  !             if (j.ne.78) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
  !            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,&
  !              1.0d0/vbld(j+nres) !d
-!            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+!            write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
              chi2=chi(itypj,itypi)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
              chi2=chi(itypj,itypi)
@@ -1530,7 +1621,7 @@
              if (rij_shift.le.0.0D0) then
                evdw=1.0D20
  !d              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
              if (rij_shift.le.0.0D0) then
                evdw=1.0D20
  !d              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-!d     &        restyp(itypi),i,restyp(itypj),j,
+!d     &        restyp(itypi,1),i,restyp(itypj,1),j,
  !d     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
                return
              endif
  !d     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
                return
              endif
@@ -1553,7 +1644,7 @@
              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
              epsi=bb**2/aa!(itypi,itypj)
              write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
              epsi=bb**2/aa!(itypi,itypj)
              write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
-              restyp(itypi),i,restyp(itypj),j, &
+              restyp(itypi,1),i,restyp(itypj,1),j, &
                epsi,sigm,chi1,chi2,chip1,chip2, &
                eps1,eps2rt**2,eps3rt**2,sig,sig0ij, &
                om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, &
                epsi,sigm,chi1,chi2,chip1,chip2, &
                eps1,eps2rt**2,eps3rt**2,sig,sig0ij, &
                om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, &
@@ -1565,6 +1656,7 @@
  !C             print *,i,j,c(1,i),c(1,j),c(2,i),c(2,j),c(3,i),c(3,j)
  !            if (energy_dec) write (iout,*) &
  !                             'evdw',i,j,evdwij
  !C             print *,i,j,c(1,i),c(1,j),c(2,i),c(2,j),c(3,i),c(3,j)
  !            if (energy_dec) write (iout,*) &
  !                             'evdw',i,j,evdwij
+!                       print *,"ZALAMKA", evdw
  
  ! Calculate gradient components.
              e1=e1*eps1*eps2rt**2*eps3rt**2
  
  ! Calculate gradient components.
              e1=e1*eps1*eps2rt**2*eps3rt**2
@@ -1597,6 +1689,7 @@
            enddo      ! j
          enddo        ! iint
        enddo          ! i
            enddo      ! j
          enddo        ! iint
        enddo          ! i
+!       print *,"ZALAMKA", evdw
  !      write (iout,*) "Number of loop steps in EGB:",ind
  !ccc      energy_dec=.false.
        return
  !      write (iout,*) "Number of loop steps in EGB:",ind
  !ccc      energy_dec=.false.
        return
@@ -1635,9 +1728,9 @@
  !     if (icall.eq.0) lprn=.true.
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !     if (icall.eq.0) lprn=.true.
  !el      ind=0
        do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
+        itypi=iabs(itype(i,1))
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=iabs(itype(i+1))
+        itypi1=iabs(itype(i+1,1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -1652,7 +1745,7 @@
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
-            itypj=iabs(itype(j))
+            itypj=iabs(itype(j,1))
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
@@ -1714,7 +1807,7 @@
              bb_aq(itypi,itypj))**(1.0D0/6.0D0)
              epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
              write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
              bb_aq(itypi,itypj))**(1.0D0/6.0D0)
              epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
              write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
-              restyp(itypi),i,restyp(itypj),j,&
+              restyp(itypi,1),i,restyp(itypj,1),j,&
                epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
                chi1,chi2,chip1,chip2,&
                eps1,eps2rt**2,eps3rt**2,&
                epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
                chi1,chi2,chip1,chip2,&
                eps1,eps2rt**2,eps3rt**2,&
@@ -1767,9 +1860,9 @@
  
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
  
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
+        itypi=iabs(itype(i,1))
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=iabs(itype(i+1))
+        itypi1=iabs(itype(i+1,1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -1780,7 +1873,7 @@
  !d        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  !d   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
  !d        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  !d   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
-            itypj=iabs(itype(j))
+            itypj=iabs(itype(j,1))
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
@@ -1854,7 +1947,7 @@
        eello_turn4=0.0d0
  !el      ind=0
        do i=iatel_s,iatel_e
        eello_turn4=0.0d0
  !el      ind=0
        do i=iatel_s,iatel_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@ -1864,7 +1957,7 @@
          num_conti=0
  !        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
          do j=ielstart(i),ielend(i)
          num_conti=0
  !        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 (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
  !el          ind=ind+1
            iteli=itel(i)
            itelj=itel(j)
  !el          ind=ind+1
            iteli=itel(i)
            itelj=itel(j)
@@ -2297,17 +2390,17 @@
          endif
  !        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
          if (i.gt. nnt+2 .and. i.lt.nct+2) then
          endif
  !        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
          if (i.gt. nnt+2 .and. i.lt.nct+2) then
-          iti = itortyp(itype(i-2))
+          iti = itortyp(itype(i-2,1))
          else
            iti=ntortyp+1
          endif
  !        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
          if (i.gt. nnt+1 .and. i.lt.nct+1) then
          else
            iti=ntortyp+1
          endif
  !        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
          if (i.gt. nnt+1 .and. i.lt.nct+1) then
-          iti1 = itortyp(itype(i-1))
+          iti1 = itortyp(itype(i-1,1))
          else
            iti1=ntortyp+1
          endif
          else
            iti1=ntortyp+1
          endif
-!          print *,iti,i,"iti",iti1,itype(i-1),itype(i-2)
+!          print *,iti,i,"iti",iti1,itype(i-1,1),itype(i-2,1)
  !d        write (iout,*) '*******i',i,' iti1',iti
  !d        write (iout,*) 'b1',b1(:,iti)
  !d        write (iout,*) 'b2',b2(:,iti)
  !d        write (iout,*) '*******i',i,' iti1',iti
  !d        write (iout,*) 'b1',b1(:,iti)
  !d        write (iout,*) 'b2',b2(:,iti)
@@ -2344,8 +2437,8 @@
          enddo
  !        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
          if (i.gt. nnt+1 .and. i.lt.nct+1) then
          enddo
  !        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
          if (i.gt. nnt+1 .and. i.lt.nct+1) then
-          if (itype(i-1).le.ntyp) then
-            iti1 = itortyp(itype(i-1))
+          if (itype(i-1,1).le.ntyp) then
+            iti1 = itortyp(itype(i-1,1))
            else
              iti1=ntortyp+1
            endif
            else
              iti1=ntortyp+1
            endif
@@ -2644,7 +2737,7 @@
  #endif
  #endif
  !d      do i=1,nres
  #endif
  #endif
  !d      do i=1,nres
-!d        iti = itortyp(itype(i))
+!d        iti = itortyp(itype(i,1))
  !d        write (iout,*) i
  !d        do j=1,2
  !d        write (iout,'(2f10.5,5x,2f10.5,5x,2f10.5)') 
  !d        write (iout,*) i
  !d        do j=1,2
  !d        write (iout,'(2f10.5,5x,2f10.5,5x,2f10.5)') 
@@ -2752,8 +2845,8 @@
  !          write (iout,*) 'i',i,' fac',fac
          enddo
        endif
  !          write (iout,*) 'i',i,' fac',fac
          enddo
        endif
-      print *,wel_loc,"wel_loc",wcorr4,wcorr5,wcorr6,wturn3,wturn4,  &
-        wturn6
+!      print *,wel_loc,"wel_loc",wcorr4,wcorr5,wcorr6,wturn3,wturn4,  &
+!        wturn6
        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
        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
@@ -2808,8 +2901,8 @@
  
  !        print *,"before iturn3 loop"
        do i=iturn3_start,iturn3_end
  
  !        print *,"before iturn3 loop"
        do i=iturn3_start,iturn3_end
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
-        .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+        .or. itype(i+2,1).eq.ntyp1 .or. itype(i+3,1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@ -2854,9 +2947,9 @@
          num_cont_hb(i)=num_conti
        enddo
        do i=iturn4_start,iturn4_end
          num_cont_hb(i)=num_conti
        enddo
        do i=iturn4_start,iturn4_end
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
-          .or. itype(i+3).eq.ntyp1 &
-          .or. itype(i+4).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+          .or. itype(i+3,1).eq.ntyp1 &
+          .or. itype(i+4,1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@ -2897,7 +2990,7 @@
  
          num_conti=num_cont_hb(i)
          call eelecij(i,i+3,ees,evdw1,eel_loc)
  
          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) &
+        if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) &
           call eturn4(i,eello_turn4)
          num_cont_hb(i)=num_conti
        enddo   ! i
           call eturn4(i,eello_turn4)
          num_cont_hb(i)=num_conti
        enddo   ! i
@@ -2905,7 +2998,7 @@
  ! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
  !
        do i=iatel_s,iatel_e
  ! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
  !
        do i=iatel_s,iatel_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@ -2947,8 +3040,8 @@
  !        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
          num_conti=num_cont_hb(i)
          do j=ielstart(i),ielend(i)
  !        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
          num_conti=num_cont_hb(i)
          do j=ielstart(i),ielend(i)
-!          write (iout,*) i,j,itype(i),itype(j)
-          if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
+!          write (iout,*) i,j,itype(i,1),itype(j,1)
+          if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
            call eelecij(i,j,ees,evdw1,eel_loc)
          enddo ! j
          num_cont_hb(i)=num_conti
            call eelecij(i,j,ees,evdw1,eel_loc)
          enddo ! j
          num_cont_hb(i)=num_conti
@@ -3450,7 +3543,7 @@
            a32=a32*fac
            a33=a33*fac
  !d          write (iout,'(4i5,4f10.5)')
            a32=a32*fac
            a33=a33*fac
  !d          write (iout,'(4i5,4f10.5)')
-!d     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
+!d     &     i,itortyp(itype(i,1)),j,itortyp(itype(j,1)),a22,a23,a32,a33
  !d          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
  !d          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
  !d     &      uy(:,j),uz(:,j)
  !d          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
  !d          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
  !d     &      uy(:,j),uz(:,j)
@@ -4283,9 +4376,9 @@
          a_temp(1,2)=a23
          a_temp(2,1)=a32
          a_temp(2,2)=a33
          a_temp(1,2)=a23
          a_temp(2,1)=a32
          a_temp(2,2)=a33
-        iti1=itortyp(itype(i+1))
-        iti2=itortyp(itype(i+2))
-        iti3=itortyp(itype(i+3))
+        iti1=itortyp(itype(i+1,1))
+        iti2=itortyp(itype(i+2,1))
+        iti3=itortyp(itype(i+3,1))
  !        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
          call transpose2(EUg(1,1,i+1),e1t(1,1))
          call transpose2(Eug(1,1,i+2),e2t(1,1))
  !        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
          call transpose2(EUg(1,1,i+1),e1t(1,1))
          call transpose2(Eug(1,1,i+2),e2t(1,1))
@@ -4554,7 +4647,7 @@
  !d    print '(a)','Enter ESCP'
  !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
  !d    print '(a)','Enter ESCP'
  !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,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))
          iteli=itel(i)
          xi=0.5D0*(c(1,i)+c(1,i+1))
          yi=0.5D0*(c(2,i)+c(2,i+1))
@@ -4563,8 +4656,8 @@
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
-          if (itype(j).eq.ntyp1) cycle
-          itypj=iabs(itype(j))
+          if (itype(j,1).eq.ntyp1) cycle
+          itypj=iabs(itype(j,1))
  ! Uncomment following three lines for SC-p interactions
  !         xj=c(1,nres+j)-xi
  !         yj=c(2,nres+j)-yi
  ! Uncomment following three lines for SC-p interactions
  !         xj=c(1,nres+j)-xi
  !         yj=c(2,nres+j)-yi
@@ -4658,7 +4751,7 @@
  !d    print '(a)','Enter ESCP'
  !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
  !d    print '(a)','Enter ESCP'
  !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,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))
          iteli=itel(i)
          xi=0.5D0*(c(1,i)+c(1,i+1))
          yi=0.5D0*(c(2,i)+c(2,i+1))
@@ -4673,7 +4766,7 @@
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
-          itypj=iabs(itype(j))
+          itypj=iabs(itype(j,1))
            if (itypj.eq.ntyp1) cycle
  ! Uncomment following three lines for SC-p interactions
  !         xj=c(1,nres+j)-xi
            if (itypj.eq.ntyp1) cycle
  ! Uncomment following three lines for SC-p interactions
  !         xj=c(1,nres+j)-xi
@@ -4849,51 +4942,107 @@
  ! 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
          if (.not.dyn_ss .and. i.le.nss) then
  ! 15/02/13 CC dynamic SSbond - additional check
  ! 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
          if (.not.dyn_ss .and. i.le.nss) then
  ! 15/02/13 CC dynamic SSbond - additional check
-         if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. &
-        iabs(itype(jjj)).eq.1) then
+         if (ii.gt.nres .and. iabs(itype(iii,1)).eq.1 .and. &
+        iabs(itype(jjj,1)).eq.1) then
            call ssbond_ene(iii,jjj,eij)
            ehpb=ehpb+2*eij
  !d          write (iout,*) "eij",eij
           endif
            call ssbond_ene(iii,jjj,eij)
            ehpb=ehpb+2*eij
  !d          write (iout,*) "eij",eij
           endif
+        else if (ii.gt.nres .and. jj.gt.nres) then
+!c Restraints from contact prediction
+          dd=dist(ii,jj)
+          if (constr_dist.eq.11) then
+            ehpb=ehpb+fordepth(i)**4.0d0 &
+               *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+            fac=fordepth(i)**4.0d0 &
+               *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+          if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, &
+            ehpb,fordepth(i),dd
+           else
+          if (dhpb1(i).gt.0.0d0) then
+            ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+            fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+!c            write (iout,*) "beta nmr",
+!c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          else
+            dd=dist(ii,jj)
+            rdis=dd-dhpb(i)
+!C Get the force constant corresponding to this distance.
+            waga=forcon(i)
+!C Calculate the contribution to energy.
+            ehpb=ehpb+waga*rdis*rdis
+!c            write (iout,*) "beta reg",dd,waga*rdis*rdis
+!C
+!C Evaluate gradient.
+!C
+            fac=waga*rdis/dd
+          endif
+          endif
+          do j=1,3
+            ggg(j)=fac*(c(j,jj)-c(j,ii))
+          enddo
+          do j=1,3
+            ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+            ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+          enddo
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
          else
          else
-! Calculate the distance between the two points and its difference from the
-! target distance.
-        dd=dist(ii,jj)
-        rdis=dd-dhpb(i)
-! Get the force constant corresponding to this distance.
-        waga=forcon(i)
-! Calculate the contribution to energy.
-        ehpb=ehpb+waga*rdis*rdis
-!
-! Evaluate gradient.
-!
-        fac=waga*rdis/dd
-!d      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
-!d   &   ' waga=',waga,' fac=',fac
-        do j=1,3
-          ggg(j)=fac*(c(j,jj)-c(j,ii))
-        enddo
-!d      print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
-! If this is a SC-SC distance, we need to calculate the contributions to the
-! Cartesian gradient in the SC vectors (ghpbx).
-        if (iii.lt.ii) then
+          dd=dist(ii,jj)
+          if (constr_dist.eq.11) then
+            ehpb=ehpb+fordepth(i)**4.0d0 &
+                *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+            fac=fordepth(i)**4.0d0 &
+                *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+          if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, &
+         ehpb,fordepth(i),dd
+           else
+          if (dhpb1(i).gt.0.0d0) then
+            ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+            fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+!c            write (iout,*) "alph nmr",
+!c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          else
+            rdis=dd-dhpb(i)
+!C Get the force constant corresponding to this distance.
+            waga=forcon(i)
+!C Calculate the contribution to energy.
+            ehpb=ehpb+waga*rdis*rdis
+!c            write (iout,*) "alpha reg",dd,waga*rdis*rdis
+!C
+!C Evaluate gradient.
+!C
+            fac=waga*rdis/dd
+          endif
+          endif
+
+            do j=1,3
+              ggg(j)=fac*(c(j,jj)-c(j,ii))
+            enddo
+!cd      print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
+!C If this is a SC-SC distance, we need to calculate the contributions to the
+!C Cartesian gradient in the SC vectors (ghpbx).
+          if (iii.lt.ii) then
            do j=1,3
              ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
              ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
            enddo
            do j=1,3
              ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
              ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
            enddo
-        endif
-!grad        do j=iii,jjj-1
-!grad          do k=1,3
-!grad            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
-!grad          enddo
-!grad        enddo
-        do k=1,3
-          ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
-          ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
-        enddo
+          endif
+!cgrad        do j=iii,jjj-1
+!cgrad          do k=1,3
+!cgrad            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
+!cgrad          enddo
+!cgrad        enddo
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
          endif
        enddo
          endif
        enddo
-      ehpb=0.5D0*ehpb
+      if (constr_dist.ne.11) ehpb=0.5D0*ehpb
+
        return
        end subroutine edis
  !-----------------------------------------------------------------------------
        return
        end subroutine edis
  !-----------------------------------------------------------------------------
@@ -4922,7 +5071,7 @@
                     deltat1,deltat2,deltat12,ed,pom1,pom2,eom1,eom2,eom12,&
                     cosphi,ggk
  
                     deltat1,deltat2,deltat12,ed,pom1,pom2,eom1,eom2,eom12,&
                     cosphi,ggk
  
-      itypi=iabs(itype(i))
+      itypi=iabs(itype(i,1))
        xi=c(1,nres+i)
        yi=c(2,nres+i)
        zi=c(3,nres+i)
        xi=c(1,nres+i)
        yi=c(2,nres+i)
        zi=c(3,nres+i)
@@ -4931,7 +5080,7 @@
        dzi=dc_norm(3,nres+i)
  !      dsci_inv=dsc_inv(itypi)
        dsci_inv=vbld_inv(nres+i)
        dzi=dc_norm(3,nres+i)
  !      dsci_inv=dsc_inv(itypi)
        dsci_inv=vbld_inv(nres+i)
-      itypj=iabs(itype(j))
+      itypj=iabs(itype(j,1))
  !      dscj_inv=dsc_inv(itypj)
        dscj_inv=vbld_inv(nres+j)
        xj=c(1,nres+j)-xi
  !      dscj_inv=dsc_inv(itypj)
        dscj_inv=vbld_inv(nres+j)
        xj=c(1,nres+j)-xi
@@ -5021,8 +5170,8 @@
  !      if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres)
  
        do i=ibondp_start,ibondp_end
  !      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
+        if (itype(i-1,1).eq.ntyp1 .and. itype(i,1).eq.ntyp1) cycle
+        if (itype(i-1,1).eq.ntyp1 .or. itype(i,1).eq.ntyp1) then
  !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          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
  !C          do j=1,3
  !C          gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
@@ -5044,11 +5193,13 @@
  !        endif
        enddo
        estr=0.5d0*AKP*estr+estr1
  !        endif
        enddo
        estr=0.5d0*AKP*estr+estr1
+!      print *,"estr_bb",estr,AKP
  !
  ! 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
  !
        do i=ibond_start,ibond_end
  !
  ! 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
  !
        do i=ibond_start,ibond_end
-        iti=iabs(itype(i))
+        iti=iabs(itype(i,1))
+        if (iti.eq.0) print *,"WARNING WRONG SETTTING",i
          if (iti.ne.10 .and. iti.ne.ntyp1) then
            nbi=nbondterm(iti)
            if (nbi.eq.1) then
          if (iti.ne.10 .and. iti.ne.ntyp1) then
            nbi=nbondterm(iti)
            if (nbi.eq.1) then
@@ -5057,6 +5208,7 @@
              "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
              AKSC(1,iti),AKSC(1,iti)*diff*diff
              estr=estr+0.5d0*AKSC(1,iti)*diff*diff
              "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
              AKSC(1,iti),AKSC(1,iti)*diff*diff
              estr=estr+0.5d0*AKSC(1,iti)*diff*diff
+!            print *,"estr_sc",estr
              do j=1,3
                gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
              enddo
              do j=1,3
                gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
              enddo
@@ -5085,6 +5237,11 @@
                usumsqder=usumsqder+ud(j)*uprod2   
              enddo
              estr=estr+uprod/usum
                usumsqder=usumsqder+ud(j)*uprod2   
              enddo
              estr=estr+uprod/usum
+!            print *,"estr_sc",estr,i
+
+             if (energy_dec) write (iout,*) &
+            "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
+            AKSC(1,iti),uprod/usum
              do j=1,3
               gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
              enddo
              do j=1,3
               gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
              enddo
@@ -5134,24 +5291,24 @@
        etheta=0.0D0
  !     write (*,'(a,i2)') 'EBEND ICG=',icg
        do i=ithet_start,ithet_end
        etheta=0.0D0
  !     write (*,'(a,i2)') 'EBEND ICG=',icg
        do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
+        if (itype(i-1,1).eq.ntyp1) cycle
  ! Zero the energy function and its derivative at 0 or pi.
          call splinthet(theta(i),0.5d0*delta,ss,ssd)
  ! Zero the energy function and its derivative at 0 or pi.
          call splinthet(theta(i),0.5d0*delta,ss,ssd)
-        it=itype(i-1)
-        ichir1=isign(1,itype(i-2))
-        ichir2=isign(1,itype(i))
-         if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1))
-         if (itype(i).eq.10) ichir2=isign(1,itype(i-1))
-         if (itype(i-1).eq.10) then
-          itype1=isign(10,itype(i-2))
-          ichir11=isign(1,itype(i-2))
-          ichir12=isign(1,itype(i-2))
-          itype2=isign(10,itype(i))
-          ichir21=isign(1,itype(i))
-          ichir22=isign(1,itype(i))
+        it=itype(i-1,1)
+        ichir1=isign(1,itype(i-2,1))
+        ichir2=isign(1,itype(i,1))
+         if (itype(i-2,1).eq.10) ichir1=isign(1,itype(i-1,1))
+         if (itype(i,1).eq.10) ichir2=isign(1,itype(i-1,1))
+         if (itype(i-1,1).eq.10) then
+          itype1=isign(10,itype(i-2,1))
+          ichir11=isign(1,itype(i-2,1))
+          ichir12=isign(1,itype(i-2,1))
+          itype2=isign(10,itype(i,1))
+          ichir21=isign(1,itype(i,1))
+          ichir22=isign(1,itype(i,1))
           endif
  
           endif
  
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+        if (i.gt.3 .and. itype(i-2,1).ne.ntyp1) then
  #ifdef OSF
            phii=phi(i)
            if (phii.ne.phii) phii=150.0
  #ifdef OSF
            phii=phi(i)
            if (phii.ne.phii) phii=150.0
@@ -5164,7 +5321,7 @@
            y(1)=0.0D0
            y(2)=0.0D0
          endif
            y(1)=0.0D0
            y(2)=0.0D0
          endif
-        if (i.lt.nres .and. itype(i).ne.ntyp1) then
+        if (i.lt.nres .and. itype(i,1).ne.ntyp1) then
  #ifdef OSF
            phii1=phi(i+1)
            if (phii1.ne.phii1) phii1=150.0
  #ifdef OSF
            phii1=phi(i+1)
            if (phii1.ne.phii1) phii1=150.0
@@ -5328,7 +5485,7 @@
        end subroutine theteng
  #else
  !-----------------------------------------------------------------------------
        end subroutine theteng
  #else
  !-----------------------------------------------------------------------------
-      subroutine ebend(etheta)
+      subroutine ebend(etheta,ethetacnstr)
  !
  ! Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
  ! angles gamma and its derivatives in consecutive thetas and gammas.
  !
  ! Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
  ! angles gamma and its derivatives in consecutive thetas and gammas.
@@ -5354,31 +5511,34 @@
  !el local variables
        integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m
        real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai
  !el local variables
        integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m
        real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai
-      real(kind=8) :: aux,etheta,ccl,ssl,scl,csl
+      real(kind=8) :: aux,etheta,ccl,ssl,scl,csl,ethetacnstr
+! local variables for constrains
+      real(kind=8) :: difi,thetiii
+       integer itheta
  
        etheta=0.0D0
        do i=ithet_start,ithet_end
  
        etheta=0.0D0
        do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
-        if (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
+        if (itype(i-1,1).eq.ntyp1) cycle
+        if (itype(i-2,1).eq.ntyp1.or.itype(i,1).eq.ntyp1) cycle
+        if (iabs(itype(i+1,1)).eq.20) iblock=2
+        if (iabs(itype(i+1,1)).ne.20) iblock=1
          dethetai=0.0d0
          dephii=0.0d0
          dephii1=0.0d0
          theti2=0.5d0*theta(i)
          dethetai=0.0d0
          dephii=0.0d0
          dephii1=0.0d0
          theti2=0.5d0*theta(i)
-        ityp2=ithetyp((itype(i-1)))
+        ityp2=ithetyp((itype(i-1,1)))
          do k=1,nntheterm
            coskt(k)=dcos(k*theti2)
            sinkt(k)=dsin(k*theti2)
          enddo
          do k=1,nntheterm
            coskt(k)=dcos(k*theti2)
            sinkt(k)=dsin(k*theti2)
          enddo
-        if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
+        if (i.gt.3 .and. itype(max0(i-3,1),1).ne.ntyp1) then
  #ifdef OSF
            phii=phi(i)
            if (phii.ne.phii) phii=150.0
  #else
            phii=phi(i)
  #endif
  #ifdef OSF
            phii=phi(i)
            if (phii.ne.phii) phii=150.0
  #else
            phii=phi(i)
  #endif
-          ityp1=ithetyp((itype(i-2)))
+          ityp1=ithetyp((itype(i-2,1)))
  ! propagation of chirality for glycine type
            do k=1,nsingle
              cosph1(k)=dcos(k*phii)
  ! propagation of chirality for glycine type
            do k=1,nsingle
              cosph1(k)=dcos(k*phii)
@@ -5386,13 +5546,13 @@
            enddo
          else
            phii=0.0d0
            enddo
          else
            phii=0.0d0
-          ityp1=ithetyp(itype(i-2))
+          ityp1=ithetyp(itype(i-2,1))
            do k=1,nsingle
              cosph1(k)=0.0d0
              sinph1(k)=0.0d0
            enddo 
          endif
            do k=1,nsingle
              cosph1(k)=0.0d0
              sinph1(k)=0.0d0
            enddo 
          endif
-        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
+        if (i.lt.nres .and. itype(i+1,1).ne.ntyp1) then
  #ifdef OSF
            phii1=phi(i+1)
            if (phii1.ne.phii1) phii1=150.0
  #ifdef OSF
            phii1=phi(i+1)
            if (phii1.ne.phii1) phii1=150.0
@@ -5400,14 +5560,14 @@
  #else
            phii1=phi(i+1)
  #endif
  #else
            phii1=phi(i+1)
  #endif
-          ityp3=ithetyp((itype(i)))
+          ityp3=ithetyp((itype(i,1)))
            do k=1,nsingle
              cosph2(k)=dcos(k*phii1)
              sinph2(k)=dsin(k*phii1)
            enddo
          else
            phii1=0.0d0
            do k=1,nsingle
              cosph2(k)=dcos(k*phii1)
              sinph2(k)=dsin(k*phii1)
            enddo
          else
            phii1=0.0d0
-          ityp3=ithetyp(itype(i))
+          ityp3=ithetyp(itype(i,1))
            do k=1,nsingle
              cosph2(k)=0.0d0
              sinph2(k)=0.0d0
            do k=1,nsingle
              cosph2(k)=0.0d0
              sinph2(k)=0.0d0
@@ -5528,6 +5688,37 @@
          if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
          gloc(nphi+i-2,icg)=wang*dethetai
        enddo
          if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
          gloc(nphi+i-2,icg)=wang*dethetai
        enddo
+!-----------thete constrains
+!      if (tor_mode.ne.2) then
+      ethetacnstr=0.0d0
+!C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+      do i=ithetaconstr_start,ithetaconstr_end
+        itheta=itheta_constr(i)
+        thetiii=theta(itheta)
+        difi=pinorm(thetiii-theta_constr0(i))
+        if (difi.gt.theta_drange(i)) then
+          difi=difi-theta_drange(i)
+          ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+          gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
+         +for_thet_constr(i)*difi**3
+        else if (difi.lt.-drange(i)) then
+          difi=difi+drange(i)
+          ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+          gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
+         +for_thet_constr(i)*difi**3
+        else
+          difi=0.0
+        endif
+       if (energy_dec) then
+        write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", &
+         i,itheta,rad2deg*thetiii, &
+         rad2deg*theta_constr0(i),  rad2deg*theta_drange(i), &
+         rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, &
+         gloc(itheta+nphi-2,icg)
+        endif
+      enddo
+!      endif
+
        return
        end subroutine ebend
  #endif
        return
        end subroutine ebend
  #endif
@@ -5565,7 +5756,7 @@
        escloc=0.0D0
  !     write (iout,'(a)') 'ESC'
        do i=loc_start,loc_end
        escloc=0.0D0
  !     write (iout,'(a)') 'ESC'
        do i=loc_start,loc_end
-        it=itype(i)
+        it=itype(i,1)
          if (it.eq.ntyp1) cycle
          if (it.eq.10) goto 1
          nlobit=nlob(iabs(it))
          if (it.eq.ntyp1) cycle
          if (it.eq.10) goto 1
          nlobit=nlob(iabs(it))
@@ -5895,7 +6086,7 @@
        delta=0.02d0*pi
        escloc=0.0D0
        do i=loc_start,loc_end
        delta=0.02d0*pi
        escloc=0.0D0
        do i=loc_start,loc_end
-        if (itype(i).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1) cycle
          costtab(i+1) =dcos(theta(i+1))
          sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
          cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
          costtab(i+1) =dcos(theta(i+1))
          sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
          cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
@@ -5904,7 +6095,7 @@
          cosfac=dsqrt(cosfac2)
          sinfac2=0.5d0/(1.0d0-costtab(i+1))
          sinfac=dsqrt(sinfac2)
          cosfac=dsqrt(cosfac2)
          sinfac2=0.5d0/(1.0d0-costtab(i+1))
          sinfac=dsqrt(sinfac2)
-        it=iabs(itype(i))
+        it=iabs(itype(i,1))
          if (it.eq.10) goto 1
  !
  !  Compute the axes of tghe local cartesian coordinates system; store in
          if (it.eq.10) goto 1
  !
  !  Compute the axes of tghe local cartesian coordinates system; store in
@@ -5922,7 +6113,7 @@
            y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
          enddo
          do j = 1,3
            y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
          enddo
          do j = 1,3
-          z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
+          z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i,1)))
          enddo     
  !       write (2,*) "i",i
  !       write (2,*) "x_prime",(x_prime(j),j=1,3)
          enddo     
  !       write (2,*) "i",i
  !       write (2,*) "x_prime",(x_prime(j),j=1,3)
@@ -5954,7 +6145,7 @@
  ! Compute the energy of the ith side cbain
  !
  !        write (2,*) "xx",xx," yy",yy," zz",zz
  ! Compute the energy of the ith side cbain
  !
  !        write (2,*) "xx",xx," yy",yy," zz",zz
-        it=iabs(itype(i))
+        it=iabs(itype(i,1))
          do j = 1,65
            x(j) = sc_parmin(j,it) 
          enddo
          do j = 1,65
            x(j) = sc_parmin(j,it) 
          enddo
@@ -5962,7 +6153,7 @@
  !c diagnostics - remove later
          xx1 = dcos(alph(2))
          yy1 = dsin(alph(2))*dcos(omeg(2))
  !c diagnostics - remove later
          xx1 = dcos(alph(2))
          yy1 = dsin(alph(2))*dcos(omeg(2))
-        zz1 = -dsign(1.0,dfloat(itype(i)))*dsin(alph(2))*dsin(omeg(2))
+        zz1 = -dsign(1.0,dfloat(itype(i,1)))*dsin(alph(2))*dsin(omeg(2))
          write(2,'(3f8.1,3f9.3,1x,3f9.3)') &
            alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,&
            xx1,yy1,zz1
          write(2,'(3f8.1,3f9.3,1x,3f9.3)') &
            alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,&
            xx1,yy1,zz1
@@ -6004,7 +6195,7 @@
  !     &   dscp1,dscp2,sumene
  !        sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
          escloc = escloc + sumene
  !     &   dscp1,dscp2,sumene
  !        sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
          escloc = escloc + sumene
-!        write (2,*) "i",i," escloc",sumene,escloc,it,itype(i)
+!        write (2,*) "i",i," escloc",sumene,escloc,it,itype(i,1)
  !     & ,zz,xx,yy
  !#define DEBUG
  #ifdef DEBUG
  !     & ,zz,xx,yy
  !#define DEBUG
  #ifdef DEBUG
@@ -6050,7 +6241,7 @@
  !        
  ! Compute the gradient of esc
  !
  !        
  ! Compute the gradient of esc
  !
-!        zz=zz*dsign(1.0,dfloat(itype(i)))
+!        zz=zz*dsign(1.0,dfloat(itype(i,1)))
          pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2
          pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2
          pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2
          pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2
          pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2
          pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2
@@ -6075,7 +6266,7 @@
                +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6) &
                +(pom1+pom2)*pom_dx
  #ifdef DEBUG
                +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6) &
                +(pom1+pom2)*pom_dx
  #ifdef DEBUG
-        write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i)
+        write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i,1)
  #endif
  !
          sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz
  #endif
  !
          sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz
@@ -6090,7 +6281,7 @@
                +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6) &
                +(pom1-pom2)*pom_dy
  #ifdef DEBUG
                +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6) &
                +(pom1-pom2)*pom_dy
  #ifdef DEBUG
-        write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i)
+        write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i,1)
  #endif
  !
          de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy &
  #endif
  !
          de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy &
@@ -6102,14 +6293,14 @@
          +x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6) &
          + ( x(14) + 2*x(17)*zz+  x(18)*xx + x(20)*yy)*(s2+s2_6)
  #ifdef DEBUG
          +x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6) &
          + ( x(14) + 2*x(17)*zz+  x(18)*xx + x(20)*yy)*(s2+s2_6)
  #ifdef DEBUG
-        write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i)
+        write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i,1)
  #endif
  !
          de_dt =  0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6) &
          -0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6) &
          +pom1*pom_dt1+pom2*pom_dt2
  #ifdef DEBUG
  #endif
  !
          de_dt =  0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6) &
          -0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6) &
          +pom1*pom_dt1+pom2*pom_dt2
  #ifdef DEBUG
-        write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i)
+        write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i,1)
  #endif
  ! 
  !
  #endif
  ! 
  !
@@ -6136,9 +6327,9 @@
           dZZ_Ci(k)=0.0d0
           do j=1,3
             dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) &
           dZZ_Ci(k)=0.0d0
           do j=1,3
             dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) &
-           *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+           *dsign(1.0d0,dfloat(itype(i,1)))*dC_norm(j,i+nres)
             dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1) &
             dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1) &
-           *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+           *dsign(1.0d0,dfloat(itype(i,1)))*dC_norm(j,i+nres)
           enddo
            
           dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
           enddo
            
           dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
@@ -6338,10 +6529,10 @@
        etors=0.0D0
        do i=iphi_start,iphi_end
        etors_ii=0.0D0
        etors=0.0D0
        do i=iphi_start,iphi_end
        etors_ii=0.0D0
-        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 &
-            .or. itype(i).eq.ntyp1) cycle
-        itori=itortyp(itype(i-2))
-        itori1=itortyp(itype(i-1))
+        if (itype(i-2,1).eq.ntyp1.or. itype(i-1,1).eq.ntyp1 &
+            .or. itype(i,1).eq.ntyp1) cycle
+        itori=itortyp(itype(i-2,1))
+        itori1=itortyp(itype(i-1,1))
          phii=phi(i)
          gloci=0.0D0
  ! Proline-Proline pair is a special case...
          phii=phi(i)
          gloci=0.0D0
  ! Proline-Proline pair is a special case...
@@ -6381,7 +6572,7 @@
               'etor',i,etors_ii
          if (lprn) &
          write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
               'etor',i,etors_ii
          if (lprn) &
          write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
-        restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,&
+        restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,itori,itori1,&
          (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
          gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
  !       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
          (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
          gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
  !       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
@@ -6441,17 +6632,17 @@
  !     lprn=.true.
        etors=0.0D0
        do i=iphi_start,iphi_end
  !     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-3).eq.ntyp1 &
-             .or. itype(i).eq.ntyp1) cycle
+        if (itype(i-2,1).eq.ntyp1 .or. itype(i-1,1).eq.ntyp1 &
+             .or. itype(i-3,1).eq.ntyp1 &
+             .or. itype(i,1).eq.ntyp1) cycle
          etors_ii=0.0D0
          etors_ii=0.0D0
-         if (iabs(itype(i)).eq.20) then
+         if (iabs(itype(i,1)).eq.20) then
           iblock=2
           else
           iblock=1
           endif
           iblock=2
           else
           iblock=1
           endif
-        itori=itortyp(itype(i-2))
-        itori1=itortyp(itype(i-1))
+        itori=itortyp(itype(i-2,1))
+        itori1=itortyp(itype(i-1,1))
          phii=phi(i)
          gloci=0.0D0
  ! Regular cosine and sine terms
          phii=phi(i)
          gloci=0.0D0
  ! Regular cosine and sine terms
@@ -6490,7 +6681,7 @@
                 'etor',i,etors_ii-v0(itori,itori1,iblock)
          if (lprn) &
          write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
                 'etor',i,etors_ii-v0(itori,itori1,iblock)
          if (lprn) &
          write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
-        restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,&
+        restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,itori,itori1,&
          (v1(j,itori,itori1,iblock),j=1,6),&
          (v2(j,itori,itori1,iblock),j=1,6)
          gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
          (v1(j,itori,itori1,iblock),j=1,6),&
          (v2(j,itori,itori1,iblock),j=1,6)
          gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
@@ -6552,18 +6743,18 @@
  !      write(iout,*) "a tu??"
        do i=iphid_start,iphid_end
          etors_d_ii=0.0D0
  !      write(iout,*) "a tu??"
        do i=iphid_start,iphid_end
          etors_d_ii=0.0D0
-        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
-            .or. itype(i-3).eq.ntyp1 &
-            .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
-        itori=itortyp(itype(i-2))
-        itori1=itortyp(itype(i-1))
-        itori2=itortyp(itype(i))
+        if (itype(i-2,1).eq.ntyp1 .or. itype(i-1,1).eq.ntyp1 &
+            .or. itype(i-3,1).eq.ntyp1 &
+            .or. itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+        itori=itortyp(itype(i-2,1))
+        itori1=itortyp(itype(i-1,1))
+        itori2=itortyp(itype(i,1))
          phii=phi(i)
          phii1=phi(i+1)
          gloci1=0.0D0
          gloci2=0.0D0
          iblock=1
          phii=phi(i)
          phii1=phi(i+1)
          gloci1=0.0D0
          gloci2=0.0D0
          iblock=1
-        if (iabs(itype(i+1)).eq.20) iblock=2
+        if (iabs(itype(i+1,1)).eq.20) iblock=2
  
  ! Regular cosine and sine terms
          do j=1,ntermd_1(itori,itori1,itori2,iblock)
  
  ! Regular cosine and sine terms
          do j=1,ntermd_1(itori,itori1,itori2,iblock)
@@ -6643,10 +6834,10 @@
  !      write (iout,*) "EBACK_SC_COR",itau_start,itau_end
        esccor=0.0D0
        do i=itau_start,itau_end
  !      write (iout,*) "EBACK_SC_COR",itau_start,itau_end
        esccor=0.0D0
        do i=itau_start,itau_end
-        if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
+        if ((itype(i-2,1).eq.ntyp1).or.(itype(i-1,1).eq.ntyp1)) cycle
          esccor_ii=0.0D0
          esccor_ii=0.0D0
-        isccori=isccortyp(itype(i-2))
-        isccori1=isccortyp(itype(i-1))
+        isccori=isccortyp(itype(i-2,1))
+        isccori1=isccortyp(itype(i-1,1))
  
  !      write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
          phii=phi(i)
  
  !      write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
          phii=phi(i)
@@ -6658,17 +6849,17 @@
  !   2 = Ca...Ca...Ca...SC
  !   3 = SC...Ca...Ca...SCi
          gloci=0.0D0
  !   2 = Ca...Ca...Ca...SC
  !   3 = SC...Ca...Ca...SCi
          gloci=0.0D0
-        if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. &
-            (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. &
-            (itype(i-1).eq.ntyp1))) &
-          .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) &
-           .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1) &
-           .or.(itype(i).eq.ntyp1))) &
-          .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. &
-            (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. &
-            (itype(i-3).eq.ntyp1)))) cycle
-        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle
-        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) &
+        if (((intertyp.eq.3).and.((itype(i-2,1).eq.10).or. &
+            (itype(i-1,1).eq.10).or.(itype(i-2,1).eq.ntyp1).or. &
+            (itype(i-1,1).eq.ntyp1))) &
+          .or. ((intertyp.eq.1).and.((itype(i-2,1).eq.10) &
+           .or.(itype(i-2,1).eq.ntyp1).or.(itype(i-1,1).eq.ntyp1) &
+           .or.(itype(i,1).eq.ntyp1))) &
+          .or.((intertyp.eq.2).and.((itype(i-1,1).eq.10).or. &
+            (itype(i-1,1).eq.ntyp1).or.(itype(i-2,1).eq.ntyp1).or. &
+            (itype(i-3,1).eq.ntyp1)))) cycle
+        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1,1).eq.ntyp1)) cycle
+        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres,1).eq.ntyp1)) &
         cycle
         do j=1,nterm_sccor(isccori,isccori1)
            v1ij=v1sccor(j,intertyp,isccori,isccori1)
         cycle
         do j=1,nterm_sccor(isccori,isccori1)
            v1ij=v1sccor(j,intertyp,isccori,isccori1)
@@ -6685,7 +6876,7 @@
          gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
          if (lprn) &
          write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
          gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
          if (lprn) &
          write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
-        restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1,&
+        restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,isccori,isccori1,&
          (v1sccor(j,intertyp,isccori,isccori1),j=1,6),&
          (v2sccor(j,intertyp,isccori,isccori1),j=1,6)
          gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
          (v1sccor(j,intertyp,isccori,isccori1),j=1,6),&
          (v2sccor(j,intertyp,isccori,isccori1),j=1,6)
          gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
@@ -7805,9 +7996,9 @@
        allocate(dipderx(3,5,4,maxconts,nres))
  !
  
        allocate(dipderx(3,5,4,maxconts,nres))
  !
  
-      iti1 = itortyp(itype(i+1))
+      iti1 = itortyp(itype(i+1,1))
        if (j.lt.nres-1) then
        if (j.lt.nres-1) then
-        itj1 = itortyp(itype(j+1))
+        itj1 = itortyp(itype(j+1,1))
        else
          itj1=ntortyp+1
        endif
        else
          itj1=ntortyp+1
        endif
@@ -7900,14 +8091,14 @@
        if (l.eq.j+1) then
  ! parallel orientation of the two CA-CA-CA frames.
          if (i.gt.1) then
        if (l.eq.j+1) then
  ! parallel orientation of the two CA-CA-CA frames.
          if (i.gt.1) then
-          iti=itortyp(itype(i))
+          iti=itortyp(itype(i,1))
          else
            iti=ntortyp+1
          endif
          else
            iti=ntortyp+1
          endif
-        itk1=itortyp(itype(k+1))
-        itj=itortyp(itype(j))
+        itk1=itortyp(itype(k+1,1))
+        itj=itortyp(itype(j,1))
          if (l.lt.nres-1) then
          if (l.lt.nres-1) then
-          itl1=itortyp(itype(l+1))
+          itl1=itortyp(itype(l+1,1))
          else
            itl1=ntortyp+1
          endif
          else
            itl1=ntortyp+1
          endif
@@ -8053,15 +8244,15 @@
        else
  ! Antiparallel orientation of the two CA-CA-CA frames.
          if (i.gt.1) then
        else
  ! Antiparallel orientation of the two CA-CA-CA frames.
          if (i.gt.1) then
-          iti=itortyp(itype(i))
+          iti=itortyp(itype(i,1))
          else
            iti=ntortyp+1
          endif
          else
            iti=ntortyp+1
          endif
-        itk1=itortyp(itype(k+1))
-        itl=itortyp(itype(l))
-        itj=itortyp(itype(j))
+        itk1=itortyp(itype(k+1,1))
+        itl=itortyp(itype(l,1))
+        itj=itortyp(itype(j,1))
          if (j.lt.nres-1) then
          if (j.lt.nres-1) then
-          itj1=itortyp(itype(j+1))
+          itj1=itortyp(itype(j+1,1))
          else 
            itj1=ntortyp+1
          endif
          else 
            itj1=ntortyp+1
          endif
@@ -8409,9 +8600,9 @@
  !d      write (iout,*)
  !d     &   'EELLO5: Contacts have occurred for peptide groups',i,j,
  !d     &   ' and',k,l
  !d      write (iout,*)
  !d     &   'EELLO5: Contacts have occurred for peptide groups',i,j,
  !d     &   ' and',k,l
-      itk=itortyp(itype(k))
-      itl=itortyp(itype(l))
-      itj=itortyp(itype(j))
+      itk=itortyp(itype(k,1))
+      itl=itortyp(itype(l,1))
+      itj=itortyp(itype(j,1))
        eello5_1=0.0d0
        eello5_2=0.0d0
        eello5_3=0.0d0
        eello5_1=0.0d0
        eello5_2=0.0d0
        eello5_3=0.0d0
@@ -8939,7 +9130,7 @@
  !       i             i                                                        C
  !                                                                              C
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !       i             i                                                        C
  !                                                                              C
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-      itk=itortyp(itype(k))
+      itk=itortyp(itype(k,1))
        s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
        s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k))
        s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k))
        s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
        s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k))
        s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k))
@@ -9235,16 +9426,16 @@
  !
  ! 4/7/01 AL Component s1 was removed, because it pertains to the respective 
  !           energy moment and not to the cluster cumulant.
  !
  ! 4/7/01 AL Component s1 was removed, because it pertains to the respective 
  !           energy moment and not to the cluster cumulant.
-      iti=itortyp(itype(i))
+      iti=itortyp(itype(i,1))
        if (j.lt.nres-1) then
        if (j.lt.nres-1) then
-        itj1=itortyp(itype(j+1))
+        itj1=itortyp(itype(j+1,1))
        else
          itj1=ntortyp+1
        endif
        else
          itj1=ntortyp+1
        endif
-      itk=itortyp(itype(k))
-      itk1=itortyp(itype(k+1))
+      itk=itortyp(itype(k,1))
+      itk1=itortyp(itype(k+1,1))
        if (l.lt.nres-1) then
        if (l.lt.nres-1) then
-        itl1=itortyp(itype(l+1))
+        itl1=itortyp(itype(l+1,1))
        else
          itl1=ntortyp+1
        endif
        else
          itl1=ntortyp+1
        endif
@@ -9356,22 +9547,22 @@
  ! 4/7/01 AL Component s1 was removed, because it pertains to the respective 
  !           energy moment and not to the cluster cumulant.
  !d      write (2,*) 'eello_graph4: wturn6',wturn6
  ! 4/7/01 AL Component s1 was removed, because it pertains to the respective 
  !           energy moment and not to the cluster cumulant.
  !d      write (2,*) 'eello_graph4: wturn6',wturn6
-      iti=itortyp(itype(i))
-      itj=itortyp(itype(j))
+      iti=itortyp(itype(i,1))
+      itj=itortyp(itype(j,1))
        if (j.lt.nres-1) then
        if (j.lt.nres-1) then
-        itj1=itortyp(itype(j+1))
+        itj1=itortyp(itype(j+1,1))
        else
          itj1=ntortyp+1
        endif
        else
          itj1=ntortyp+1
        endif
-      itk=itortyp(itype(k))
+      itk=itortyp(itype(k,1))
        if (k.lt.nres-1) then
        if (k.lt.nres-1) then
-        itk1=itortyp(itype(k+1))
+        itk1=itortyp(itype(k+1,1))
        else
          itk1=ntortyp+1
        endif
        else
          itk1=ntortyp+1
        endif
-      itl=itortyp(itype(l))
+      itl=itortyp(itype(l,1))
        if (l.lt.nres-1) then
        if (l.lt.nres-1) then
-        itl1=itortyp(itype(l+1))
+        itl1=itortyp(itype(l+1,1))
        else
          itl1=ntortyp+1
        endif
        else
          itl1=ntortyp+1
        endif
@@ -9598,11 +9789,11 @@
        j=i+4
        k=i+1
        l=i+3
        j=i+4
        k=i+1
        l=i+3
-      iti=itortyp(itype(i))
-      itk=itortyp(itype(k))
-      itk1=itortyp(itype(k+1))
-      itl=itortyp(itype(l))
-      itj=itortyp(itype(j))
+      iti=itortyp(itype(i,1))
+      itk=itortyp(itype(k,1))
+      itk1=itortyp(itype(k+1,1))
+      itl=itortyp(itype(l,1))
+      itj=itortyp(itype(j,1))
  !d      write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj
  !d      write (2,*) 'i',i,' k',k,' j',j,' l',l
  !d      if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then
  !d      write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj
  !d      write (2,*) 'i',i,' k',k,' j',j,' l',l
  !d      if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then
@@ -10141,15 +10332,14 @@
                        wturn6*gcorr6_turn_long(j,i)+ &
                        wstrain*ghpbc(j,i) &
                       +wliptran*gliptranc(j,i) &
                        wturn6*gcorr6_turn_long(j,i)+ &
                        wstrain*ghpbc(j,i) &
                       +wliptran*gliptranc(j,i) &
+                     +gradafm(j,i) &
                       +welec*gshieldc(j,i) &
                       +wcorr*gshieldc_ec(j,i) &
                       +wturn3*gshieldc_t3(j,i)&
                       +wturn4*gshieldc_t4(j,i)&
                       +wel_loc*gshieldc_ll(j,i)&
                       +welec*gshieldc(j,i) &
                       +wcorr*gshieldc_ec(j,i) &
                       +wturn3*gshieldc_t3(j,i)&
                       +wturn4*gshieldc_t4(j,i)&
                       +wel_loc*gshieldc_ll(j,i)&
-                     +wtube*gg_tube(j,i)
- 
-
-
+                     +wtube*gg_tube(j,i) &
+                     +wbond_nucl*gradb_nucl(j,i)
          enddo
        enddo 
  #else
          enddo
        enddo 
  #else
@@ -10166,13 +10356,13 @@
                        wturn6*gcorr6_turn_long(j,i)+ &
                        wstrain*ghpbc(j,i) &
                       +wliptran*gliptranc(j,i) &
                        wturn6*gcorr6_turn_long(j,i)+ &
                        wstrain*ghpbc(j,i) &
                       +wliptran*gliptranc(j,i) &
+                     +gradafm(j,i) &
                       +welec*gshieldc(j,i)&
                       +wcorr*gshieldc_ec(j,i) &
                       +wturn4*gshieldc_t4(j,i) &
                       +wel_loc*gshieldc_ll(j,i)&
                       +welec*gshieldc(j,i)&
                       +wcorr*gshieldc_ec(j,i) &
                       +wturn4*gshieldc_t4(j,i) &
                       +wel_loc*gshieldc_ll(j,i)&
-                     +wtube*gg_tube(j,i)
-
-
+                     +wtube*gg_tube(j,i) &
+                     +wbond_nucl*gradb_nucl(j,i)
  
          enddo
        enddo 
  
          enddo
        enddo 
@@ -10318,6 +10508,7 @@
                        wsccor*gsccorc(j,i) &
                       +wscloc*gscloc(j,i)  &
                       +wliptran*gliptranc(j,i) &
                        wsccor*gsccorc(j,i) &
                       +wscloc*gscloc(j,i)  &
                       +wliptran*gliptranc(j,i) &
+                     +gradafm(j,i) &
                       +welec*gshieldc(j,i) &
                       +welec*gshieldc_loc(j,i) &
                       +wcorr*gshieldc_ec(j,i) &
                       +welec*gshieldc(j,i) &
                       +welec*gshieldc_loc(j,i) &
                       +wcorr*gshieldc_ec(j,i) &
@@ -10328,7 +10519,9 @@
                       +wturn4*gshieldc_loc_t4(j,i) &
                       +wel_loc*gshieldc_ll(j,i) &
                       +wel_loc*gshieldc_loc_ll(j,i) &
                       +wturn4*gshieldc_loc_t4(j,i) &
                       +wel_loc*gshieldc_ll(j,i) &
                       +wel_loc*gshieldc_loc_ll(j,i) &
-                     +wtube*gg_tube(j,i)
+                     +wtube*gg_tube(j,i) &
+                     +wbond_nucl*gradb_nucl(j,i)
+
  
  
  #else
  
  
  #else
@@ -10350,6 +10543,7 @@
                        wturn6*gcorr6_turn(j,i)+ &
                        wsccor*gsccorc(j,i) &
                       +wscloc*gscloc(j,i) &
                        wturn6*gcorr6_turn(j,i)+ &
                        wsccor*gsccorc(j,i) &
                       +wscloc*gscloc(j,i) &
+                     +gradafm(j,i) &
                       +wliptran*gliptranc(j,i) &
                       +welec*gshieldc(j,i) &
                       +welec*gshieldc_loc(j,) &
                       +wliptran*gliptranc(j,i) &
                       +welec*gshieldc(j,i) &
                       +welec*gshieldc_loc(j,) &
@@ -10361,7 +10555,9 @@
                       +wturn4*gshieldc_loc_t4(j,i) &
                       +wel_loc*gshieldc_ll(j,i) &
                       +wel_loc*gshieldc_loc_ll(j,i) &
                       +wturn4*gshieldc_loc_t4(j,i) &
                       +wel_loc*gshieldc_ll(j,i) &
                       +wel_loc*gshieldc_loc_ll(j,i) &
-                     +wtube*gg_tube(j,i)
+                     +wtube*gg_tube(j,i) &
+                     +wbond_nucl*gradb_nucl(j,i) 
+
  
  
  
  
  
  
@@ -10377,7 +10573,9 @@
                         +wturn3*gshieldx_t3(j,i) &
                         +wturn4*gshieldx_t4(j,i) &
                         +wel_loc*gshieldx_ll(j,i)&
                         +wturn3*gshieldx_t3(j,i) &
                         +wturn4*gshieldx_t4(j,i) &
                         +wel_loc*gshieldx_ll(j,i)&
-                       +wtube*gg_tube_sc(j,i)
+                       +wtube*gg_tube_sc(j,i)   &
+                       +wbond_nucl*gradbx_nucl(j,i) 
+
  
  
          enddo
  
  
          enddo
@@ -10587,7 +10785,7 @@
  !      include 'COMMON.CALC'
  !      include 'COMMON.IOUNITS'
        real(kind=8), dimension(3) :: dcosom1,dcosom2
  !      include 'COMMON.CALC'
  !      include 'COMMON.IOUNITS'
        real(kind=8), dimension(3) :: dcosom1,dcosom2
-
+!      print *,"wchodze"
        eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
        eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
        eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
        eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
        eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
        eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
@@ -10963,7 +11161,7 @@
  ! Derivatives in alpha and omega:
  !
        do i=2,nres-1
  ! Derivatives in alpha and omega:
  !
        do i=2,nres-1
-!       dsci=dsc(itype(i))
+!       dsci=dsc(itype(i,1))
          dsci=vbld(i+nres)
  #ifdef OSF
          alphi=alph(i)
          dsci=vbld(i+nres)
  #ifdef OSF
          alphi=alph(i)
@@ -12051,9 +12249,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
  !      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -12064,7 +12262,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !d        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  !d   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
  !d        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  !d   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
@@ -12141,9 +12339,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
  !      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -12156,7 +12354,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !d        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  !d   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
  !d        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  !d   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
@@ -12231,9 +12429,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
  !     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -12242,7 +12440,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
@@ -12262,7 +12460,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !d            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
  !d            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
  !d            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
  !d            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
  !d            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
  !d            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d   &          restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d   &          restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
  !d   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
  !d   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
  !d   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
  !d   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
  !d   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
  !d   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
@@ -12318,9 +12516,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
  !     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -12329,7 +12527,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
@@ -12349,7 +12547,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !d            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
  !d            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
  !d            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
  !d            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
  !d            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
  !d            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d   &          restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d   &          restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
  !d   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
  !d   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
  !d   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
  !d   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
  !d   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
  !d   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
@@ -12417,9 +12615,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     endif
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !     endif
  !el      ind=0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -12434,7 +12632,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
@@ -12475,7 +12673,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
  !d              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
  !d              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
-!d     &          restyp(itypi),i,restyp(itypj),j,
+!d     &          restyp(itypi,1),i,restyp(itypj,1),j,
  !d     &          epsi,sigm,chi1,chi2,chip1,chip2,
  !d     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
  !d     &          om1,om2,om12,1.0D0/dsqrt(rrij),
  !d     &          epsi,sigm,chi1,chi2,chip1,chip2,
  !d     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
  !d     &          om1,om2,om12,1.0D0/dsqrt(rrij),
@@ -12537,9 +12735,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     endif
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !     endif
  !el      ind=0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -12554,7 +12752,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
@@ -12595,7 +12793,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
  !d              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
  !d              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
-!d     &          restyp(itypi),i,restyp(itypj),j,
+!d     &          restyp(itypi,1),i,restyp(itypj,1),j,
  !d     &          epsi,sigm,chi1,chi2,chip1,chip2,
  !d     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
  !d     &          om1,om2,om12,1.0D0/dsqrt(rrij),
  !d     &          epsi,sigm,chi1,chi2,chip1,chip2,
  !d     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
  !d     &          om1,om2,om12,1.0D0/dsqrt(rrij),
@@ -12657,9 +12855,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     if (icall.eq.0) lprn=.false.
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !     if (icall.eq.0) lprn=.false.
  !el      ind=0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -12705,21 +12903,43 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
              IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
              IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
-              call dyn_ssbond_ene(i,j,evdwij)
-              evdw=evdw+evdwij
-              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
-                              'evdw',i,j,evdwij,' ss'
+!              call dyn_ssbond_ene(i,j,evdwij)
+!              evdw=evdw+evdwij
+!              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+!                              'evdw',i,j,evdwij,' ss'
  !              if (energy_dec) write (iout,*) &
  !                              'evdw',i,j,evdwij,' ss'
  !              if (energy_dec) write (iout,*) &
  !                              'evdw',i,j,evdwij,' ss'
+!             do k=j+1,iend(i,iint)
+!C search over all next residues
+!              if (dyn_ss_mask(k)) then
+!C check if they are cysteins
+!C              write(iout,*) 'k=',k
+
+!c              write(iout,*) "PRZED TRI", evdwij
+!               evdwij_przed_tri=evdwij
+!              call triple_ssbond_ene(i,j,k,evdwij)
+!c               if(evdwij_przed_tri.ne.evdwij) then
+!c                 write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+!c               endif
+
+!c              write(iout,*) "PO TRI", evdwij
+!C call the energy function that removes the artifical triple disulfide
+!C bond the soubroutine is located in ssMD.F
+!              evdw=evdw+evdwij
+              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+                            'evdw',i,j,evdwij,'tss'
+!              endif!dyn_ss_mask(k)
+!             enddo! k
+
              ELSE
  !el            ind=ind+1
              ELSE
  !el            ind=ind+1
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
  !            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
  !     &       1.0d0/vbld(j+nres)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
  !            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
  !     &       1.0d0/vbld(j+nres)
-!            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+!            write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
              chi2=chi(itypj,itypi)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
              chi2=chi(itypj,itypi)
@@ -12823,7 +13043,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                if (rij_shift.le.0.0D0) then
                  evdw=1.0D20
  !d                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
                if (rij_shift.le.0.0D0) then
                  evdw=1.0D20
  !d                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-!d     &          restyp(itypi),i,restyp(itypj),j,
+!d     &          restyp(itypi,1),i,restyp(itypj,1),j,
  !d     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
                  return
                endif
  !d     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
                  return
                endif
@@ -12844,7 +13064,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
                write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
                write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
-                restyp(itypi),i,restyp(itypj),j,&
+                restyp(itypi,1),i,restyp(itypj,1),j,&
                  epsi,sigm,chi1,chi2,chip1,chip2,&
                  eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
                  om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
                  epsi,sigm,chi1,chi2,chip1,chip2,&
                  eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
                  om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
@@ -12915,9 +13135,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     if (icall.eq.0) lprn=.false.
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !     if (icall.eq.0) lprn=.false.
  !el      ind=0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -12973,17 +13193,39 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                evdw=evdw+evdwij
                if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
                                'evdw',i,j,evdwij,' ss'
                evdw=evdw+evdwij
                if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
                                'evdw',i,j,evdwij,' ss'
+             do k=j+1,iend(i,iint)
+!C search over all next residues
+              if (dyn_ss_mask(k)) then
+!C check if they are cysteins
+!C              write(iout,*) 'k=',k
+
+!c              write(iout,*) "PRZED TRI", evdwij
+!               evdwij_przed_tri=evdwij
+              call triple_ssbond_ene(i,j,k,evdwij)
+!c               if(evdwij_przed_tri.ne.evdwij) then
+!c                 write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+!c               endif
+
+!c              write(iout,*) "PO TRI", evdwij
+!C call the energy function that removes the artifical triple disulfide
+!C bond the soubroutine is located in ssMD.F
+              evdw=evdw+evdwij
+              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+                            'evdw',i,j,evdwij,'tss'
+              endif!dyn_ss_mask(k)
+             enddo! k
+
  !              if (energy_dec) write (iout,*) &
  !                              'evdw',i,j,evdwij,' ss'
              ELSE
  !el            ind=ind+1
  !              if (energy_dec) write (iout,*) &
  !                              'evdw',i,j,evdwij,' ss'
              ELSE
  !el            ind=ind+1
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
  !            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
  !     &       1.0d0/vbld(j+nres)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
  !            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
  !     &       1.0d0/vbld(j+nres)
-!            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+!            write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
              chi2=chi(itypj,itypi)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
              chi2=chi(itypj,itypi)
@@ -13092,7 +13334,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                if (rij_shift.le.0.0D0) then
                  evdw=1.0D20
  !d                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
                if (rij_shift.le.0.0D0) then
                  evdw=1.0D20
  !d                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-!d     &          restyp(itypi),i,restyp(itypj),j,
+!d     &          restyp(itypi,1),i,restyp(itypj,1),j,
  !d     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
                  return
                endif
  !d     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
                  return
                endif
@@ -13113,7 +13355,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
                write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
                write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
-                restyp(itypi),i,restyp(itypj),j,&
+                restyp(itypi,1),i,restyp(itypj,1),j,&
                  epsi,sigm,chi1,chi2,chip1,chip2,&
                  eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
                  om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
                  epsi,sigm,chi1,chi2,chip1,chip2,&
                  eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
                  om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
@@ -13183,9 +13425,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     if (icall.eq.0) lprn=.true.
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !     if (icall.eq.0) lprn=.true.
  !el      ind=0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -13200,7 +13442,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
@@ -13256,7 +13498,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
                write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
                write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
-                restyp(itypi),i,restyp(itypj),j,&
+                restyp(itypi,1),i,restyp(itypj,1),j,&
                  epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
                  chi1,chi2,chip1,chip2,&
                  eps1,eps2rt**2,eps3rt**2,&
                  epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
                  chi1,chi2,chip1,chip2,&
                  eps1,eps2rt**2,eps3rt**2,&
@@ -13312,9 +13554,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     if (icall.eq.0) lprn=.true.
  !el      ind=0
        do i=iatsc_s,iatsc_e
  !     if (icall.eq.0) lprn=.true.
  !el      ind=0
        do i=iatsc_s,iatsc_e
-        itypi=itype(i)
+        itypi=itype(i,1)
          if (itypi.eq.ntyp1) cycle
          if (itypi.eq.ntyp1) cycle
-        itypi1=itype(i+1)
+        itypi1=itype(i+1,1)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@ -13329,7 +13571,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
  !el            ind=ind+1
-            itypj=itype(j)
+            itypj=itype(j,1)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              if (itypj.eq.ntyp1) cycle
  !            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
@@ -13385,7 +13627,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
                write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
                sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
                epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
                write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
-                restyp(itypi),i,restyp(itypj),j,&
+                restyp(itypi,1),i,restyp(itypj,1),j,&
                  epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
                  chi1,chi2,chip1,chip2,&
                  eps1,eps2rt**2,eps3rt**2,&
                  epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
                  chi1,chi2,chip1,chip2,&
                  eps1,eps2rt**2,eps3rt**2,&
@@ -13536,8 +13778,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  ! Loop over i,i+2 and i,i+3 pairs of the peptide groups
  !
        do i=iturn3_start,iturn3_end
  ! Loop over i,i+2 and i,i+3 pairs of the peptide groups
  !
        do i=iturn3_start,iturn3_end
-        if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 &
-        .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1.or. itype(i+1,1).eq.ntyp1 &
+        .or. itype(i+2,1).eq.ntyp1 .or. itype(i+3,1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@ -13559,9 +13801,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          num_cont_hb(i)=num_conti
        enddo
        do i=iturn4_start,iturn4_end
          num_cont_hb(i)=num_conti
        enddo
        do i=iturn4_start,iturn4_end
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
-          .or. itype(i+3).eq.ntyp1 &
-          .or. itype(i+4).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+          .or. itype(i+3,1).eq.ntyp1 &
+          .or. itype(i+4,1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@ -13579,7 +13821,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
            if (zmedi.lt.0) zmedi=zmedi+boxzsize
          num_conti=num_cont_hb(i)
          call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
            if (zmedi.lt.0) zmedi=zmedi+boxzsize
          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.ntyp1) &
+        if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) &
            call eturn4(i,eello_turn4)
          num_cont_hb(i)=num_conti
        enddo   ! i
            call eturn4(i,eello_turn4)
          num_cont_hb(i)=num_conti
        enddo   ! i
@@ -13587,7 +13829,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  ! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
  !
        do i=iatel_s,iatel_e
  ! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
  !
        do i=iatel_s,iatel_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@ -13606,7 +13848,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
          num_conti=num_cont_hb(i)
          do j=ielstart(i),ielend(i)
  !        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.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+          if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
            call eelecij_scale(i,j,ees,evdw1,eel_loc)
          enddo ! j
          num_cont_hb(i)=num_conti
            call eelecij_scale(i,j,ees,evdw1,eel_loc)
          enddo ! j
          num_cont_hb(i)=num_conti
@@ -13990,7 +14232,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
            a32=a32*fac
            a33=a33*fac
  !d          write (iout,'(4i5,4f10.5)')
            a32=a32*fac
            a33=a33*fac
  !d          write (iout,'(4i5,4f10.5)')
-!d     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
+!d     &     i,itortyp(itype(i,1)),j,itortyp(itype(j,1)),a22,a23,a32,a33
  !d          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
  !d          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
  !d     &      uy(:,j),uz(:,j)
  !d          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
  !d          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
  !d     &      uy(:,j),uz(:,j)
@@ -14467,7 +14709,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     & " iatel_e_vdw",iatel_e_vdw
        call flush(iout)
        do i=iatel_s_vdw,iatel_e_vdw
  !     & " iatel_e_vdw",iatel_e_vdw
        call flush(iout)
        do i=iatel_s_vdw,iatel_e_vdw
-        if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1.or. itype(i+1,1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@ -14488,7 +14730,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !     &   ' ielend',ielend_vdw(i)
          call flush(iout)
          do j=ielstart_vdw(i),ielend_vdw(i)
  !     &   ' ielend',ielend_vdw(i)
          call flush(iout)
          do j=ielstart_vdw(i),ielend_vdw(i)
-          if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+          if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
  !el          ind=ind+1
            iteli=itel(i)
            itelj=itel(j)
  !el          ind=ind+1
            iteli=itel(i)
            itelj=itel(j)
@@ -14621,7 +14863,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !d    print '(a)','Enter ESCP'
  !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
  !d    print '(a)','Enter ESCP'
  !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,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))
          iteli=itel(i)
          xi=0.5D0*(c(1,i)+c(1,i+1))
          yi=0.5D0*(c(2,i)+c(2,i+1))
@@ -14636,7 +14878,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
-          itypj=itype(j)
+          itypj=itype(j,1)
            if (itypj.eq.ntyp1) cycle
  ! Uncomment following three lines for SC-p interactions
  !         xj=c(1,nres+j)-xi
            if (itypj.eq.ntyp1) cycle
  ! Uncomment following three lines for SC-p interactions
  !         xj=c(1,nres+j)-xi
@@ -14780,7 +15022,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !d    print '(a)','Enter ESCP'
  !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
  !d    print '(a)','Enter ESCP'
  !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1 .or. itype(i+1,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))
          iteli=itel(i)
          xi=0.5D0*(c(1,i)+c(1,i+1))
          yi=0.5D0*(c(2,i)+c(2,i+1))
@@ -14795,7 +15037,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
-          itypj=itype(j)
+          itypj=itype(j,1)
            if (itypj.eq.ntyp1) cycle
  ! Uncomment following three lines for SC-p interactions
  !         xj=c(1,nres+j)-xi
            if (itypj.eq.ntyp1) cycle
  ! Uncomment following three lines for SC-p interactions
  !         xj=c(1,nres+j)-xi
@@ -15255,7 +15497,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !el local variables
        integer :: i,nres6
        real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,esccor,etors_d,etors
  !el local variables
        integer :: i,nres6
        real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,esccor,etors_d,etors
-      real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr
+      real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr,ethetacnstr
        nres6=6*nres
  
  !      write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
        nres6=6*nres
  
  !      write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
@@ -15410,7 +15652,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !
  ! Calculate the virtual-bond-angle energy.
  !
  !
  ! Calculate the virtual-bond-angle energy.
  !
-      call ebend(ebe)
+      call ebend(ebe,ethetacnstr)
  !
  ! Calculate the SC local energy.
  !
  !
  ! Calculate the SC local energy.
  !
@@ -15496,7 +15738,35 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        endif
        return
        end function gnmr1prim
        endif
        return
        end function gnmr1prim
-!-----------------------------------------------------------------------------
+!----------------------------------------------------------------------------
+      real(kind=8) function rlornmr1(y,ymin,ymax,sigma)
+      real(kind=8) y,ymin,ymax,sigma
+      real(kind=8) wykl /4.0d0/
+      if (y.lt.ymin) then
+        rlornmr1=(ymin-y)**wykl/((ymin-y)**wykl+sigma**wykl)
+      else if (y.gt.ymax) then
+        rlornmr1=(y-ymax)**wykl/((y-ymax)**wykl+sigma**wykl)
+      else
+        rlornmr1=0.0d0
+      endif
+      return
+      end function rlornmr1
+!------------------------------------------------------------------------------
+      real(kind=8) function rlornmr1prim(y,ymin,ymax,sigma)
+      real(kind=8) y,ymin,ymax,sigma
+      real(kind=8) wykl /4.0d0/
+      if (y.lt.ymin) then
+        rlornmr1prim=-(ymin-y)**(wykl-1)*sigma**wykl*wykl/ &
+        ((ymin-y)**wykl+sigma**wykl)**2
+      else if (y.gt.ymax) then
+        rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/ &
+        ((y-ymax)**wykl+sigma**wykl)**2
+      else
+        rlornmr1prim=0.0d0
+      endif
+      return
+      end function rlornmr1prim
+
        real(kind=8) function harmonic(y,ymax)
  !      implicit none
        real(kind=8) :: y,ymax
        real(kind=8) function harmonic(y,ymax)
  !      implicit none
        real(kind=8) :: y,ymax
@@ -15591,7 +15861,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        enddo
        if (n.le.nphi+ntheta) goto 10
        do i=2,nres-1
        enddo
        if (n.le.nphi+ntheta) goto 10
        do i=2,nres-1
-       if (itype(i).ne.10) then
+       if (itype(i,1).ne.10) then
            galphai=0.0D0
           gomegai=0.0D0
           do k=1,3
            galphai=0.0D0
           gomegai=0.0D0
           do k=1,3
@@ -15886,6 +16156,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
            gshieldc_loc_ll(j,i)=0.0d0
            gg_tube(j,i)=0.0d0
            gg_tube_sc(j,i)=0.0d0
            gshieldc_loc_ll(j,i)=0.0d0
            gg_tube(j,i)=0.0d0
            gg_tube_sc(j,i)=0.0d0
+          gradafm(j,i)=0.0d0
+          gradb_nucl(j,i)=0.0d0
+          gradbx_nucl(j,i)=0.0d0
            do intertyp=1,3
             gloc_sc(intertyp,i,icg)=0.0d0
            enddo
            do intertyp=1,3
             gloc_sc(intertyp,i,icg)=0.0d0
            enddo
@@ -16019,10 +16292,10 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          do j=1,3
            dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/&
           vbld(i-1)
          do j=1,3
            dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/&
           vbld(i-1)
-          if (itype(i-1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint
+          if (itype(i-1,1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint
            dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/&
           vbld(i)
            dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/&
           vbld(i)
-          if (itype(i-1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint
+          if (itype(i-1,1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint
          enddo
        enddo
  #if defined(MPI) && defined(PARINTDER)
          enddo
        enddo
  #if defined(MPI) && defined(PARINTDER)
@@ -16031,7 +16304,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  #else
        do i=3,nres
  #endif
  #else
        do i=3,nres
  #endif
-      if ((itype(i-1).ne.10).and.(itype(i-1).ne.ntyp1)) then
+      if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1)) then
          cost1=dcos(omicron(1,i))
          sint1=sqrt(1-cost1*cost1)
          cost2=dcos(omicron(2,i))
          cost1=dcos(omicron(1,i))
          sint1=sqrt(1-cost1*cost1)
          cost2=dcos(omicron(2,i))
@@ -16055,7 +16328,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
            dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) &
             +cost2*(-dc_norm(j,i-1+nres)))/ &
            vbld(i-1+nres)
            dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) &
             +cost2*(-dc_norm(j,i-1+nres)))/ &
            vbld(i-1+nres)
-!          write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
+!          write(iout,*) "vbld", i,itype(i,1),vbld(i-1+nres)
            domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
          enddo
         endif
            domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
          enddo
         endif
@@ -16070,7 +16343,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  #else
        do i=4,nres      
  #endif
  #else
        do i=4,nres      
  #endif
-!        if (itype(i-1).eq.21 .or. itype(i-2).eq.21 ) cycle
+!        if (itype(i-1,1).eq.21 .or. itype(i-2,1).eq.21 ) cycle
  ! the conventional case
          sint=dsin(theta(i))
         sint1=dsin(theta(i-1))
  ! the conventional case
          sint=dsin(theta(i))
         sint1=dsin(theta(i-1))
@@ -16095,7 +16368,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
              ctgt=cost/sint
              ctgt1=cost1/sint1
              cosg_inv=1.0d0/cosg
              ctgt=cost/sint
              ctgt1=cost1/sint1
              cosg_inv=1.0d0/cosg
-            if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
+            if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
             dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) &
                -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
              dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
             dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) &
                -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
              dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
@@ -16113,7 +16386,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !   Obtaining the gamma derivatives from cosine derivative
          else
             do j=1,3
  !   Obtaining the gamma derivatives from cosine derivative
          else
             do j=1,3
-           if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
+           if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
             dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* &
            dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* &
             dc_norm(j,i-3))/vbld(i-2)
             dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* &
            dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* &
             dc_norm(j,i-3))/vbld(i-2)
@@ -16137,9 +16410,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        do i=3,nres
  !elwrite(iout,*) " vecpr",i,nres
  #endif
        do i=3,nres
  !elwrite(iout,*) " vecpr",i,nres
  #endif
-       if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
-!       if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or.
-!     &     (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
+       if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle
+!       if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10).or.
+!     &     (itype(i-1,1).eq.ntyp1).or.(itype(i,1).eq.ntyp1)) cycle
  !c dtauangle(j,intertyp,dervityp,residue number)
  !c INTERTYP=1 SC...Ca...Ca..Ca
  ! the conventional case
  !c dtauangle(j,intertyp,dervityp,residue number)
  !c INTERTYP=1 SC...Ca...Ca..Ca
  ! the conventional case
@@ -16217,8 +16490,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  #else
        do i=4,nres
  #endif
  #else
        do i=4,nres
  #endif
-       if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. &
-          (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) cycle
+       if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. &
+          (itype(i-2,1).eq.ntyp1).or.(itype(i-3,1).eq.ntyp1)) cycle
  ! the conventional case
          sint=dsin(omicron(1,i))
          sint1=dsin(theta(i-1))
  ! the conventional case
          sint=dsin(omicron(1,i))
          sint1=dsin(theta(i-1))
@@ -16292,8 +16565,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        do i=3,nres
  #endif
  ! the conventional case
        do i=3,nres
  #endif
  ! the conventional case
-      if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. &
-      (itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
+      if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. &
+      (itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle
          sint=dsin(omicron(1,i))
          sint1=dsin(omicron(2,i-1))
          sing=dsin(tauangle(3,i))
          sint=dsin(omicron(1,i))
          sint1=dsin(omicron(2,i-1))
          sing=dsin(tauangle(3,i))
@@ -16363,7 +16636,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  #else
          do i=2,nres-1          
  #endif
  #else
          do i=2,nres-1          
  #endif
-          if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then        
+          if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then    
               fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
               fac6=fac5/vbld(i)
               fac7=fac5*fac5
               fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
               fac6=fac5/vbld(i)
               fac7=fac5*fac5
@@ -16597,7 +16870,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        write (iout,*) &
         "Analytical (upper) and numerical (lower) gradient of alpha"
        do i=2,nres-1
        write (iout,*) &
         "Analytical (upper) and numerical (lower) gradient of alpha"
        do i=2,nres-1
-       if(itype(i).ne.10) then
+       if(itype(i,1).ne.10) then
                     do j=1,3
               dcji=dc(j,i-1)
                       dc(j,i-1)=dcji+aincr
                     do j=1,3
               dcji=dc(j,i-1)
                       dc(j,i-1)=dcji+aincr
@@ -16633,7 +16906,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        write (iout,*) &
         "Analytical (upper) and numerical (lower) gradient of omega"
        do i=2,nres-1
        write (iout,*) &
         "Analytical (upper) and numerical (lower) gradient of omega"
        do i=2,nres-1
-       if(itype(i).ne.10) then
+       if(itype(i,1).ne.10) then
                     do j=1,3
               dcji=dc(j,i-1)
                       dc(j,i-1)=dcji+aincr
                     do j=1,3
               dcji=dc(j,i-1)
                       dc(j,i-1)=dcji+aincr
@@ -16699,7 +16972,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                         (cref(3,jl,kkk)-cref(3,il,kkk))**2)
              dij=dist(il,jl)
              qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
                         (cref(3,jl,kkk)-cref(3,il,kkk))**2)
              dij=dist(il,jl)
              qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
-            if (itype(il).ne.10 .or. itype(jl).ne.10) then
+            if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
                nl=nl+1
                d0ijCM=dsqrt( &
                       (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
                nl=nl+1
                d0ijCM=dsqrt( &
                       (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
@@ -16726,7 +16999,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                         (cref(3,jl,kkk)-cref(3,il,kkk))**2)
              dij=dist(il,jl)
              qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
                         (cref(3,jl,kkk)-cref(3,il,kkk))**2)
              dij=dist(il,jl)
              qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
-            if (itype(il).ne.10 .or. itype(jl).ne.10) then
+            if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
                nl=nl+1
                d0ijCM=dsqrt( &
                       (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
                nl=nl+1
                d0ijCM=dsqrt( &
                       (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
@@ -16788,7 +17061,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                dqwol(k,jl)=dqwol(k,jl)-ddqij
              enddo
                      
                dqwol(k,jl)=dqwol(k,jl)-ddqij
              enddo
                      
-            if (itype(il).ne.10 .or. itype(jl).ne.10) then
+            if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
                nl=nl+1
                d0ijCM=dsqrt( &
                       (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
                nl=nl+1
                d0ijCM=dsqrt( &
                       (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
@@ -16829,7 +17102,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
                dqwol(k,il)=dqwol(k,il)+ddqij
                dqwol(k,jl)=dqwol(k,jl)-ddqij
              enddo
                dqwol(k,il)=dqwol(k,il)+ddqij
                dqwol(k,jl)=dqwol(k,jl)-ddqij
              enddo
-            if (itype(il).ne.10 .or. itype(jl).ne.10) then
+            if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
                nl=nl+1
                d0ijCM=dsqrt( &
                       (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
                nl=nl+1
                d0ijCM=dsqrt( &
                       (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
@@ -17320,13 +17593,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !el      allocate(dyn_ssbond_ij(iatsc_s:iatsc_e,nres))
  !el      allocate(dyn_ssbond_ij(0:nres+4,nres))
  
  !el      allocate(dyn_ssbond_ij(iatsc_s:iatsc_e,nres))
  !el      allocate(dyn_ssbond_ij(0:nres+4,nres))
  
-      itypi=itype(i)
+      itypi=itype(i,1)
        dxi=dc_norm(1,nres+i)
        dyi=dc_norm(2,nres+i)
        dzi=dc_norm(3,nres+i)
        dsci_inv=vbld_inv(i+nres)
  
        dxi=dc_norm(1,nres+i)
        dyi=dc_norm(2,nres+i)
        dzi=dc_norm(3,nres+i)
        dsci_inv=vbld_inv(i+nres)
  
-      itypj=itype(j)
+      itypj=itype(j,1)
        xj=c(1,nres+j)-c(1,nres+i)
        yj=c(2,nres+j)-c(2,nres+i)
        zj=c(3,nres+j)-c(3,nres+i)
        xj=c(1,nres+j)-c(1,nres+i)
        yj=c(2,nres+j)-c(2,nres+i)
        zj=c(3,nres+j)-c(3,nres+i)
@@ -17637,6 +17910,176 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  
        return
        end subroutine dyn_ssbond_ene
  
        return
        end subroutine dyn_ssbond_ene
+!--------------------------------------------------------------------------
+         subroutine triple_ssbond_ene(resi,resj,resk,eij)
+!      implicit none
+!      Includes
+      use calc_data
+      use comm_sschecks
+!      include 'DIMENSIONS'
+!      include 'COMMON.SBRIDGE'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.DERIV'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.VAR'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+       use MD_data
+!      include 'COMMON.MD'
+!      use MD, only: totT,t_bath
+#endif
+#endif
+      double precision h_base
+      external h_base
+
+!c     Input arguments
+      integer resi,resj,resk,m,itypi,itypj,itypk
+
+!c     Output arguments
+      double precision eij,eij1,eij2,eij3
+
+!c     Local variables
+      logical havebond
+!c      integer itypi,itypj,k,l
+      double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
+      double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
+      double precision xik,yik,zik,xjk,yjk,zjk,dxk,dyk,dzk
+      double precision sig0ij,ljd,sig,fac,e1,e2
+      double precision dcosom1(3),dcosom2(3),ed
+      double precision pom1,pom2
+      double precision ljA,ljB,ljXs
+      double precision d_ljB(1:3)
+      double precision ssA,ssB,ssC,ssXs
+      double precision ssxm,ljxm,ssm,ljm
+      double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+      eij=0.0
+      if (dtriss.eq.0) return
+      i=resi
+      j=resj
+      k=resk
+!C      write(iout,*) resi,resj,resk
+      itypi=itype(i,1)
+      dxi=dc_norm(1,nres+i)
+      dyi=dc_norm(2,nres+i)
+      dzi=dc_norm(3,nres+i)
+      dsci_inv=vbld_inv(i+nres)
+      xi=c(1,nres+i)
+      yi=c(2,nres+i)
+      zi=c(3,nres+i)
+      itypj=itype(j,1)
+      xj=c(1,nres+j)
+      yj=c(2,nres+j)
+      zj=c(3,nres+j)
+
+      dxj=dc_norm(1,nres+j)
+      dyj=dc_norm(2,nres+j)
+      dzj=dc_norm(3,nres+j)
+      dscj_inv=vbld_inv(j+nres)
+      itypk=itype(k,1)
+      xk=c(1,nres+k)
+      yk=c(2,nres+k)
+      zk=c(3,nres+k)
+
+      dxk=dc_norm(1,nres+k)
+      dyk=dc_norm(2,nres+k)
+      dzk=dc_norm(3,nres+k)
+      dscj_inv=vbld_inv(k+nres)
+      xij=xj-xi
+      xik=xk-xi
+      xjk=xk-xj
+      yij=yj-yi
+      yik=yk-yi
+      yjk=yk-yj
+      zij=zj-zi
+      zik=zk-zi
+      zjk=zk-zj
+      rrij=(xij*xij+yij*yij+zij*zij)
+      rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
+      rrik=(xik*xik+yik*yik+zik*zik)
+      rik=dsqrt(rrik)
+      rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
+      rjk=dsqrt(rrjk)
+!C there are three combination of distances for each trisulfide bonds
+!C The first case the ith atom is the center
+!C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
+!C distance y is second distance the a,b,c,d are parameters derived for
+!C this problem d parameter was set as a penalty currenlty set to 1.
+      if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then
+      eij1=0.0d0
+      else
+      eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss)
+      endif
+!C second case jth atom is center
+      if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then
+      eij2=0.0d0
+      else
+      eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss)
+      endif
+!C the third case kth atom is the center
+      if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then
+      eij3=0.0d0
+      else
+      eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss)
+      endif
+!C      eij2=0.0
+!C      eij3=0.0
+!C      eij1=0.0
+      eij=eij1+eij2+eij3
+!C      write(iout,*)i,j,k,eij
+!C The energy penalty calculated now time for the gradient part 
+!C derivative over rij
+      fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+      -eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)
+            gg(1)=xij*fac/rij
+            gg(2)=yij*fac/rij
+            gg(3)=zij*fac/rij
+      do m=1,3
+        gvdwx(m,i)=gvdwx(m,i)-gg(m)
+        gvdwx(m,j)=gvdwx(m,j)+gg(m)
+      enddo
+
+      do l=1,3
+        gvdwc(l,i)=gvdwc(l,i)-gg(l)
+        gvdwc(l,j)=gvdwc(l,j)+gg(l)
+      enddo
+!C now derivative over rik
+      fac=-eij1**2/dtriss* &
+      (-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+      -eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+            gg(1)=xik*fac/rik
+            gg(2)=yik*fac/rik
+            gg(3)=zik*fac/rik
+      do m=1,3
+        gvdwx(m,i)=gvdwx(m,i)-gg(m)
+        gvdwx(m,k)=gvdwx(m,k)+gg(m)
+      enddo
+      do l=1,3
+        gvdwc(l,i)=gvdwc(l,i)-gg(l)
+        gvdwc(l,k)=gvdwc(l,k)+gg(l)
+      enddo
+!C now derivative over rjk
+      fac=-eij2**2/dtriss* &
+      (-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)- &
+      eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+            gg(1)=xjk*fac/rjk
+            gg(2)=yjk*fac/rjk
+            gg(3)=zjk*fac/rjk
+      do m=1,3
+        gvdwx(m,j)=gvdwx(m,j)-gg(m)
+        gvdwx(m,k)=gvdwx(m,k)+gg(m)
+      enddo
+      do l=1,3
+        gvdwc(l,j)=gvdwc(l,j)-gg(l)
+        gvdwc(l,k)=gvdwc(l,k)+gg(l)
+      enddo
+      return
+      end subroutine triple_ssbond_ene
+
+
+
  !-----------------------------------------------------------------------------
        real(kind=8) function h_base(x,deriv)
  !     A smooth function going 0->1 in range [0,1]
  !-----------------------------------------------------------------------------
        real(kind=8) function h_base(x,deriv)
  !     A smooth function going 0->1 in range [0,1]
@@ -17783,15 +18226,18 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        diff=newnss-nss
  
  !mc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
        diff=newnss-nss
  
  !mc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
-
+!       print *,newnss,nss,maxdim
        do i=1,nss
          found=.false.
        do i=1,nss
          found=.false.
+!        print *,newnss
          do j=1,newnss
          do j=1,newnss
+!!          print *,j
            if (idssb(i).eq.newihpb(j) .and. &
                 jdssb(i).eq.newjhpb(j)) found=.true.
          enddo
  #ifndef CLUST
  #ifndef WHAM
            if (idssb(i).eq.newihpb(j) .and. &
                 jdssb(i).eq.newjhpb(j)) found=.true.
          enddo
  #ifndef CLUST
  #ifndef WHAM
+!        write(iout,*) "found",found,i,j
          if (.not.found.and.fg_rank.eq.0) &
              write(iout,'(a15,f12.2,f8.1,2i5)') &
               "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
          if (.not.found.and.fg_rank.eq.0) &
              write(iout,'(a15,f12.2,f8.1,2i5)') &
               "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
@@ -17802,11 +18248,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        do i=1,newnss
          found=.false.
          do j=1,nss
        do i=1,newnss
          found=.false.
          do j=1,nss
+!          print *,i,j
            if (newihpb(i).eq.idssb(j) .and. &
                 newjhpb(i).eq.jdssb(j)) found=.true.
          enddo
  #ifndef CLUST
  #ifndef WHAM
            if (newihpb(i).eq.idssb(j) .and. &
                 newjhpb(i).eq.jdssb(j)) found=.true.
          enddo
  #ifndef CLUST
  #ifndef WHAM
+!        write(iout,*) "found",found,i,j
          if (.not.found.and.fg_rank.eq.0) &
              write(iout,'(a15,f12.2,f8.1,2i5)') &
               "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
          if (.not.found.and.fg_rank.eq.0) &
              write(iout,'(a15,f12.2,f8.1,2i5)') &
               "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
@@ -17836,10 +18284,10 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        real(kind=8) :: fracinbuf,eliptran,sslip,positi,ssgradlip
        integer :: i
        eliptran=0.0
        real(kind=8) :: fracinbuf,eliptran,sslip,positi,ssgradlip
        integer :: i
        eliptran=0.0
-      print *, "I am in eliptran"
+!      print *, "I am in eliptran"
        do i=ilip_start,ilip_end
  !C       do i=1,1
        do i=ilip_start,ilip_end
  !C       do i=1,1
-        if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1).or.(i.eq.nres))&
+        if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1).or.(i.eq.nres))&
           cycle
  
          positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
           cycle
  
          positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
@@ -17885,7 +18333,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         enddo
  ! here starts the side chain transfer
         do i=ilip_start,ilip_end
         enddo
  ! here starts the side chain transfer
         do i=ilip_start,ilip_end
-        if (itype(i).eq.ntyp1) cycle
+        if (itype(i,1).eq.ntyp1) cycle
          positi=(mod(c(3,i+nres),boxzsize))
          if (positi.le.0) positi=positi+boxzsize
  !C       print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
          positi=(mod(c(3,i+nres),boxzsize))
          if (positi.le.0) positi=positi+boxzsize
  !C       print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
@@ -17901,25 +18349,25 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !C lipbufthick is thickenes of lipid buffore
           sslip=sscalelip(fracinbuf)
           ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
  !C lipbufthick is thickenes of lipid buffore
           sslip=sscalelip(fracinbuf)
           ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
-         eliptran=eliptran+sslip*liptranene(itype(i))
+         eliptran=eliptran+sslip*liptranene(itype(i,1))
           gliptranx(3,i)=gliptranx(3,i) &
           gliptranx(3,i)=gliptranx(3,i) &
-      +ssgradlip*liptranene(itype(i))
+      +ssgradlip*liptranene(itype(i,1))
           gliptranc(3,i-1)= gliptranc(3,i-1) &
           gliptranc(3,i-1)= gliptranc(3,i-1) &
-      +ssgradlip*liptranene(itype(i))
+      +ssgradlip*liptranene(itype(i,1))
  !C         print *,"doing sccale for lower part"
          elseif (positi.gt.bufliptop) then
           fracinbuf=1.0d0-  &
        ((bordliptop-positi)/lipbufthick)
           sslip=sscalelip(fracinbuf)
           ssgradlip=sscagradlip(fracinbuf)/lipbufthick
  !C         print *,"doing sccale for lower part"
          elseif (positi.gt.bufliptop) then
           fracinbuf=1.0d0-  &
        ((bordliptop-positi)/lipbufthick)
           sslip=sscalelip(fracinbuf)
           ssgradlip=sscagradlip(fracinbuf)/lipbufthick
-         eliptran=eliptran+sslip*liptranene(itype(i))
+         eliptran=eliptran+sslip*liptranene(itype(i,1))
           gliptranx(3,i)=gliptranx(3,i)  &
           gliptranx(3,i)=gliptranx(3,i)  &
-       +ssgradlip*liptranene(itype(i))
+       +ssgradlip*liptranene(itype(i,1))
           gliptranc(3,i-1)= gliptranc(3,i-1) &
           gliptranc(3,i-1)= gliptranc(3,i-1) &
-      +ssgradlip*liptranene(itype(i))
+      +ssgradlip*liptranene(itype(i,1))
  !C          print *, "doing sscalefor top part",sslip,fracinbuf
          else
  !C          print *, "doing sscalefor top part",sslip,fracinbuf
          else
-         eliptran=eliptran+liptranene(itype(i))
+         eliptran=eliptran+liptranene(itype(i,1))
  !C         print *,"I am in true lipid"
          endif
          endif ! if in lipid or buffor
  !C         print *,"I am in true lipid"
          endif
          endif ! if in lipid or buffor
@@ -17942,7 +18390,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !C and r0 is the excluded size of nanotube (can be set to 0 if we want just a 
  !C simple Kihara potential
        subroutine calctube(Etube)
  !C and r0 is the excluded size of nanotube (can be set to 0 if we want just a 
  !C simple Kihara potential
        subroutine calctube(Etube)
-      real(kind=8) :: vectube(3),enetube(nres*2)
+      real(kind=8),dimension(3) :: vectube
        real(kind=8) :: Etube,xtemp,xminact,yminact,& 
         ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, &
         sc_aa_tube,sc_bb_tube
        real(kind=8) :: Etube,xtemp,xminact,yminact,& 
         ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, &
         sc_aa_tube,sc_bb_tube
@@ -17956,7 +18404,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !C for UNRES
         do i=itube_start,itube_end
  !C lets ommit dummy atoms for now
  !C for UNRES
         do i=itube_start,itube_end
  !C lets ommit dummy atoms for now
-       if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+       if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
  !C now calculate distance from center of tube and direction vectors
        xmin=boxxsize
        ymin=boxysize
  !C now calculate distance from center of tube and direction vectors
        xmin=boxxsize
        ymin=boxysize
@@ -18019,7 +18467,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  
         do i=itube_start,itube_end
  !C Lets not jump over memory as we use many times iti
  
         do i=itube_start,itube_end
  !C Lets not jump over memory as we use many times iti
-         iti=itype(i)
+         iti=itype(i,1)
  !C lets ommit dummy atoms for now
           if ((iti.eq.ntyp1)  &
  !C in UNRES uncomment the line below as GLY has no side-chain...
  !C lets ommit dummy atoms for now
           if ((iti.eq.ntyp1)  &
  !C in UNRES uncomment the line below as GLY has no side-chain...
@@ -18103,7 +18551,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !C and r0 is the excluded size of nanotube (can be set to 0 if we want just a 
  !C simple Kihara potential
        subroutine calctube2(Etube)
  !C and r0 is the excluded size of nanotube (can be set to 0 if we want just a 
  !C simple Kihara potential
        subroutine calctube2(Etube)
-      real(kind=8) :: vectube(3),enetube(nres*2)
+            real(kind=8),dimension(3) :: vectube
        real(kind=8) :: Etube,xtemp,xminact,yminact,&
         ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi,fracinbuf,&
         sstube,ssgradtube,sc_aa_tube,sc_bb_tube
        real(kind=8) :: Etube,xtemp,xminact,yminact,&
         ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi,fracinbuf,&
         sstube,ssgradtube,sc_aa_tube,sc_bb_tube
@@ -18119,7 +18567,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         do i=itube_start,itube_end
  !C lets ommit dummy atoms for now
  
         do i=itube_start,itube_end
  !C lets ommit dummy atoms for now
  
-       if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+       if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
  !C now calculate distance from center of tube and direction vectors
  !C      vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
  !C          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
  !C now calculate distance from center of tube and direction vectors
  !C      vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
  !C          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
@@ -18180,12 +18628,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !C lipbufthick is thickenes of lipid buffore
           sstube=sscalelip(fracinbuf)
           ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
  !C lipbufthick is thickenes of lipid buffore
           sstube=sscalelip(fracinbuf)
           ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
-!C         print *,ssgradtube, sstube,tubetranene(itype(i))
+!C         print *,ssgradtube, sstube,tubetranene(itype(i,1))
           enetube(i)=enetube(i)+sstube*tubetranenepep
  !C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
           enetube(i)=enetube(i)+sstube*tubetranenepep
  !C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C     &+ssgradtube*tubetranene(itype(i))
+!C     &+ssgradtube*tubetranene(itype(i,1))
  !C         gg_tube(3,i-1)= gg_tube(3,i-1)
  !C         gg_tube(3,i-1)= gg_tube(3,i-1)
-!C     &+ssgradtube*tubetranene(itype(i))
+!C     &+ssgradtube*tubetranene(itype(i,1))
  !C         print *,"doing sccale for lower part"
          elseif (positi.gt.buftubetop) then
           fracinbuf=1.0d0-  &
  !C         print *,"doing sccale for lower part"
          elseif (positi.gt.buftubetop) then
           fracinbuf=1.0d0-  &
@@ -18194,9 +18642,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           ssgradtube=sscagradlip(fracinbuf)/tubebufthick
           enetube(i)=enetube(i)+sstube*tubetranenepep
  !C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
           ssgradtube=sscagradlip(fracinbuf)/tubebufthick
           enetube(i)=enetube(i)+sstube*tubetranenepep
  !C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C     &+ssgradtube*tubetranene(itype(i))
+!C     &+ssgradtube*tubetranene(itype(i,1))
  !C         gg_tube(3,i-1)= gg_tube(3,i-1)
  !C         gg_tube(3,i-1)= gg_tube(3,i-1)
-!C     &+ssgradtube*tubetranene(itype(i))
+!C     &+ssgradtube*tubetranene(itype(i,1))
  !C          print *, "doing sscalefor top part",sslip,fracinbuf
          else
           sstube=1.0d0
  !C          print *, "doing sscalefor top part",sslip,fracinbuf
          else
           sstube=1.0d0
@@ -18237,7 +18685,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !C        print *,gg_tube(1,0),"TU"
          do i=itube_start,itube_end
  !C Lets not jump over memory as we use many times iti
  !C        print *,gg_tube(1,0),"TU"
          do i=itube_start,itube_end
  !C Lets not jump over memory as we use many times iti
-         iti=itype(i)
+         iti=itype(i,1)
  !C lets ommit dummy atoms for now
           if ((iti.eq.ntyp1) &
  !!C in UNRES uncomment the line below as GLY has no side-chain...
  !C lets ommit dummy atoms for now
           if ((iti.eq.ntyp1) &
  !!C in UNRES uncomment the line below as GLY has no side-chain...
@@ -18269,12 +18717,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !C lipbufthick is thickenes of lipid buffore
           sstube=sscalelip(fracinbuf)
           ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
  !C lipbufthick is thickenes of lipid buffore
           sstube=sscalelip(fracinbuf)
           ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
-!C         print *,ssgradtube, sstube,tubetranene(itype(i))
-         enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+!C         print *,ssgradtube, sstube,tubetranene(itype(i,1))
+         enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
  !C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
  !C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C     &+ssgradtube*tubetranene(itype(i))
+!C     &+ssgradtube*tubetranene(itype(i,1))
  !C         gg_tube(3,i-1)= gg_tube(3,i-1)
  !C         gg_tube(3,i-1)= gg_tube(3,i-1)
-!C     &+ssgradtube*tubetranene(itype(i))
+!C     &+ssgradtube*tubetranene(itype(i,1))
  !C         print *,"doing sccale for lower part"
          elseif (positi.gt.buftubetop) then
           fracinbuf=1.0d0- &
  !C         print *,"doing sccale for lower part"
          elseif (positi.gt.buftubetop) then
           fracinbuf=1.0d0- &
@@ -18282,16 +18730,16 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  
           sstube=sscalelip(fracinbuf)
           ssgradtube=sscagradlip(fracinbuf)/tubebufthick
  
           sstube=sscalelip(fracinbuf)
           ssgradtube=sscagradlip(fracinbuf)/tubebufthick
-         enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+         enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
  !C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
  !C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C     &+ssgradtube*tubetranene(itype(i))
+!C     &+ssgradtube*tubetranene(itype(i,1))
  !C         gg_tube(3,i-1)= gg_tube(3,i-1)
  !C         gg_tube(3,i-1)= gg_tube(3,i-1)
-!C     &+ssgradtube*tubetranene(itype(i))
+!C     &+ssgradtube*tubetranene(itype(i,1))
  !C          print *, "doing sscalefor top part",sslip,fracinbuf
          else
           sstube=1.0d0
           ssgradtube=0.0d0
  !C          print *, "doing sscalefor top part",sslip,fracinbuf
          else
           sstube=1.0d0
           ssgradtube=0.0d0
-         enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+         enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
  !C         print *,"I am in true lipid"
          endif
          else
  !C         print *,"I am in true lipid"
          endif
          else
@@ -18340,15 +18788,15 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
          end subroutine calctube2
  !=====================================================================================================================================
        subroutine calcnano(Etube)
          end subroutine calctube2
  !=====================================================================================================================================
        subroutine calcnano(Etube)
-      real(kind=8) :: vectube(3),enetube(nres*2), &
-      enecavtube(nres*2)
+      real(kind=8),dimension(3) :: vectube
+      
        real(kind=8) :: Etube,xtemp,xminact,yminact,&
         ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,denominator,faccav,&
         sc_aa_tube,sc_bb_tube,zmin,ztemp,zminact
        real(kind=8) :: Etube,xtemp,xminact,yminact,&
         ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,denominator,faccav,&
         sc_aa_tube,sc_bb_tube,zmin,ztemp,zminact
-       integer:: i,j,iti
+       integer:: i,j,iti,r
  
        Etube=0.0d0
  
        Etube=0.0d0
-      print *,itube_start,itube_end,"poczatek"
+!      print *,itube_start,itube_end,"poczatek"
        do i=itube_start,itube_end
          enetube(i)=0.0d0
          enetube(i+nres)=0.0d0
        do i=itube_start,itube_end
          enetube(i)=0.0d0
          enetube(i+nres)=0.0d0
@@ -18358,7 +18806,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !C for UNRES
         do i=itube_start,itube_end
  !C lets ommit dummy atoms for now
  !C for UNRES
         do i=itube_start,itube_end
  !C lets ommit dummy atoms for now
-       if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+       if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
  !C now calculate distance from center of tube and direction vectors
        xmin=boxxsize
        ymin=boxysize
  !C now calculate distance from center of tube and direction vectors
        xmin=boxxsize
        ymin=boxysize
@@ -18441,7 +18889,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !C         fac=fac+faccav
  !C 667     continue
           endif
  !C         fac=fac+faccav
  !C 667     continue
           endif
-
+          if (energy_dec) write(iout,*),i,rdiff,enetube(i),enecavtube(i)
          do j=1,3
          gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
          gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
          do j=1,3
          gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
          gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
@@ -18451,7 +18899,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         do i=itube_start,itube_end
          enecavtube(i)=0.0d0
  !C Lets not jump over memory as we use many times iti
         do i=itube_start,itube_end
          enecavtube(i)=0.0d0
  !C Lets not jump over memory as we use many times iti
-         iti=itype(i)
+         iti=itype(i,1)
  !C lets ommit dummy atoms for now
           if ((iti.eq.ntyp1) &
  !C in UNRES uncomment the line below as GLY has no side-chain...
  !C lets ommit dummy atoms for now
           if ((iti.eq.ntyp1) &
  !C in UNRES uncomment the line below as GLY has no side-chain...
@@ -18544,6 +18992,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
            gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
            gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
           enddo
            gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
            gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
           enddo
+          if (energy_dec) write(iout,*),i,rdiff,enetube(i+nres),enecavtube(i+nres)
          enddo
  
  
          enddo
  
  
@@ -18552,6 +19001,23 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
            Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) &
           +enecavtube(i+nres)
          enddo
            Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) &
           +enecavtube(i+nres)
          enddo
+!        do i=1,20
+!         print *,"begin", i,"a"
+!         do r=1,10000
+!          rdiff=r/100.0d0
+!          rdiff6=rdiff**6.0d0
+!          sc_aa_tube=sc_aa_tube_par(i)
+!          sc_bb_tube=sc_bb_tube_par(i)
+!          enetube(i)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+!          denominator=(1.0d0+dcavtub(i)*rdiff6*rdiff6)
+!          enecavtube(i)=   &
+!         (bcavtub(i)*rdiff+acavtub(i)*dsqrt(rdiff)+ccavtub(i)) &
+!         /denominator
+
+!          print '(5(f10.3,1x))',rdiff,enetube(i),enecavtube(i),enecavtube(i)+enetube(i)
+!         enddo
+!         print *,"end",i,"a"
+!        enddo
  !C        print *,"ETUBE", etube
          return
          end subroutine calcnano
  !C        print *,"ETUBE", etube
          return
          end subroutine calcnano
@@ -18582,14 +19048,14 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        enddo
        do i=ivec_start,ivec_end
  !C      do i=1,nres-1
        enddo
        do i=ivec_start,ivec_end
  !C      do i=1,nres-1
-!C      if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+!C      if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle
        ishield_list(i)=0
        ishield_list(i)=0
-      if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+      if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle
  !Cif there two consequtive dummy atoms there is no peptide group between them
  !C the line below has to be changed for FGPROC>1
        VolumeTotal=0.0
        do k=1,nres
  !Cif there two consequtive dummy atoms there is no peptide group between them
  !C the line below has to be changed for FGPROC>1
        VolumeTotal=0.0
        do k=1,nres
-       if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+       if ((itype(k,1).eq.ntyp1).or.(itype(k,1).eq.10)) cycle
         dist_pep_side=0.0
         dist_side_calf=0.0
         do j=1,3
         dist_pep_side=0.0
         dist_side_calf=0.0
         do j=1,3
@@ -18644,8 +19110,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           enddo
          endif
  !C this is what is now we have the distance scaling now volume...
           enddo
          endif
  !C this is what is now we have the distance scaling now volume...
-      short=short_r_sidechain(itype(k))
-      long=long_r_sidechain(itype(k))
+      short=short_r_sidechain(itype(k,1))
+      long=long_r_sidechain(itype(k,1))
        costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
        sinthet=short/dist_pep_side*costhet
  !C now costhet_grad
        costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
        sinthet=short/dist_pep_side*costhet
  !C now costhet_grad
@@ -18731,10 +19197,66 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        enddo
        fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
       
        enddo
        fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
       
-!C      write(2,*) "TOTAL VOLUME",i,itype(i),fac_shield(i)
+!C      write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i)
        enddo
        return
        end subroutine set_shield_fac2
        enddo
        return
        end subroutine set_shield_fac2
+!----------------------------------------------------------------------------
+! SOUBROUTINE FOR AFM
+       subroutine AFMvel(Eafmforce)
+       use MD_data, only:totTafm
+      real(kind=8),dimension(3) :: diffafm
+      real(kind=8) :: afmdist,Eafmforce
+       integer :: i
+!C Only for check grad COMMENT if not used for checkgrad
+!C      totT=3.0d0
+!C--------------------------------------------------------
+!C      print *,"wchodze"
+      afmdist=0.0d0
+      Eafmforce=0.0d0
+      do i=1,3
+      diffafm(i)=c(i,afmend)-c(i,afmbeg)
+      afmdist=afmdist+diffafm(i)**2
+      enddo
+      afmdist=dsqrt(afmdist)
+!      totTafm=3.0
+      Eafmforce=0.5d0*forceAFMconst &
+      *(distafminit+totTafm*velAFMconst-afmdist)**2
+!C      Eafmforce=-forceAFMconst*(dist-distafminit)
+      do i=1,3
+      gradafm(i,afmend-1)=-forceAFMconst* &
+       (distafminit+totTafm*velAFMconst-afmdist) &
+       *diffafm(i)/afmdist
+      gradafm(i,afmbeg-1)=forceAFMconst* &
+      (distafminit+totTafm*velAFMconst-afmdist) &
+      *diffafm(i)/afmdist
+      enddo
+!      print *,'AFM',Eafmforce,totTafm*velAFMconst,afmdist
+      return
+      end subroutine AFMvel
+!---------------------------------------------------------
+       subroutine AFMforce(Eafmforce)
+
+      real(kind=8),dimension(3) :: diffafm
+!      real(kind=8) ::afmdist
+      real(kind=8) :: afmdist,Eafmforce
+      integer :: i
+      afmdist=0.0d0
+      Eafmforce=0.0d0
+      do i=1,3
+      diffafm(i)=c(i,afmend)-c(i,afmbeg)
+      afmdist=afmdist+diffafm(i)**2
+      enddo
+      afmdist=dsqrt(afmdist)
+!      print *,afmdist,distafminit
+      Eafmforce=-forceAFMconst*(afmdist-distafminit)
+      do i=1,3
+      gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/afmdist
+      gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/afmdist
+      enddo
+!C      print *,'AFM',Eafmforce
+      return
+      end subroutine AFMforce
  
  !-----------------------------------------------------------------------------
  #ifdef WHAM
  
  !-----------------------------------------------------------------------------
  #ifdef WHAM
@@ -18982,6 +19504,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        allocate(grad_shield(3,-1:nres))
        allocate(gg_tube_sc(3,-1:nres))
        allocate(gg_tube(3,-1:nres))
        allocate(grad_shield(3,-1:nres))
        allocate(gg_tube_sc(3,-1:nres))
        allocate(gg_tube(3,-1:nres))
+      allocate(gradafm(3,-1:nres))
+      allocate(gradb_nucl(3,-1:nres))
+      allocate(gradbx_nucl(3,-1:nres))
  !(3,maxres)
        allocate(grad_shield_side(3,50,nres))
        allocate(grad_shield_loc(3,50,nres))
  !(3,maxres)
        allocate(grad_shield_side(3,50,nres))
        allocate(grad_shield_loc(3,50,nres))
@@ -19101,14 +19626,18 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  !        enddo
  !      enddo
  
  !        enddo
  !      enddo
  
-      if (nss.gt.0) then
-        allocate(idssb(nss),jdssb(nss))
+!      if (nss.gt.0) then
+        allocate(idssb(maxdim),jdssb(maxdim))
+!        allocate(newihpb(nss),newjhpb(nss))
  !(maxdim)
  !(maxdim)
-      endif
+!      endif
        allocate(ishield_list(nres))
        allocate(shield_list(50,nres))
        allocate(dyn_ss_mask(nres))
        allocate(fac_shield(nres))
        allocate(ishield_list(nres))
        allocate(shield_list(50,nres))
        allocate(dyn_ss_mask(nres))
        allocate(fac_shield(nres))
+      allocate(enetube(nres*2))
+      allocate(enecavtube(nres*2))
+
  !(maxres)
        dyn_ss_mask(:)=.false.
  !----------------------
  !(maxres)
        dyn_ss_mask(:)=.false.
  !----------------------
@@ -19155,6 +19684,283 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
  
        return
        end subroutine alloc_ener_arrays
  
        return
        end subroutine alloc_ener_arrays
+!-----------------------------------------------------------------
+      subroutine ebond_nucl(estr_nucl)
+!c
+!c Evaluate the energy of stretching of the CA-CA and CA-SC virtual bonds
+!c 
+      
+      real(kind=8),dimension(3) :: u,ud
+      real(kind=8) :: usum,uprod,uprod1,uprod2,usumsqder
+      real(kind=8) :: estr_nucl,diff
+      integer :: iti,i,j,k,nbi
+      estr_nucl=0.0d0
+!C      print *,"I enter ebond"
+      if (energy_dec) &
+      write (iout,*) "ibondp_start,ibondp_end",&
+       ibondp_nucl_start,ibondp_nucl_end
+      do i=ibondp_nucl_start,ibondp_nucl_end
+        if (itype(i-1,2).eq.ntyp1_molec(2) .or. &
+         itype(i,2).eq.ntyp1_molec(2)) cycle
+!          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)
+
+          diff = vbld(i)-vbldp0_nucl
+          if(energy_dec)write(iout,*) "estr_nucl_bb" , i,vbld(i),&
+          vbldp0_nucl,diff,AKP_nucl*diff*diff
+          estr_nucl=estr_nucl+diff*diff
+          print *,estr_nucl
+          do j=1,3
+            gradb_nucl(j,i-1)=AKP_nucl*diff*dc(j,i-1)/vbld(i)
+          enddo
+!c          write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
+      enddo
+      estr_nucl=0.5d0*AKP_nucl*estr_nucl
+      print *,"partial sum", estr_nucl,AKP_nucl
+
+      if (energy_dec) &
+      write (iout,*) "ibondp_start,ibondp_end",&
+       ibond_nucl_start,ibond_nucl_end
+
+      do i=ibond_nucl_start,ibond_nucl_end
+!C        print *, "I am stuck",i
+        iti=itype(i,2)
+        if (iti.eq.ntyp1_molec(2)) cycle
+          nbi=nbondterm_nucl(iti)
+!C        print *,iti,nbi
+          if (nbi.eq.1) then
+            diff=vbld(i+nres)-vbldsc0_nucl(1,iti)
+
+            if (energy_dec) &
+           write (iout,*) "estr_nucl_sc", i,iti,vbld(i+nres),vbldsc0_nucl(1,iti),diff, &
+           AKSC_nucl(1,iti),AKSC_nucl(1,iti)*diff*diff
+            estr_nucl=estr_nucl+0.5d0*AKSC_nucl(1,iti)*diff*diff
+            print *,estr_nucl
+            do j=1,3
+              gradbx_nucl(j,i)=AKSC_nucl(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
+            enddo
+          else
+            do j=1,nbi
+              diff=vbld(i+nres)-vbldsc0_nucl(j,iti)
+              ud(j)=aksc_nucl(j,iti)*diff
+              u(j)=abond0_nucl(j,iti)+0.5d0*ud(j)*diff
+            enddo
+            uprod=u(1)
+            do j=2,nbi
+              uprod=uprod*u(j)
+            enddo
+            usum=0.0d0
+            usumsqder=0.0d0
+            do j=1,nbi
+              uprod1=1.0d0
+              uprod2=1.0d0
+              do k=1,nbi
+                if (k.ne.j) then
+                  uprod1=uprod1*u(k)
+                  uprod2=uprod2*u(k)*u(k)
+                endif
+              enddo
+              usum=usum+uprod1
+              usumsqder=usumsqder+ud(j)*uprod2
+            enddo
+            estr_nucl=estr_nucl+uprod/usum
+            do j=1,3
+             gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
+            enddo
+        endif
+      enddo
+!C      print *,"I am about to leave ebond"
+      return
+      end subroutine ebond_nucl
+
+!-----------------------------------------------------------------------------
+      subroutine ebend_nucl(etheta_nucl)
+      real(kind=8),dimension(nntheterm_nucl+1) :: coskt,sinkt !mmaxtheterm
+      real(kind=8),dimension(nsingle_nucl+1) :: cosph1,sinph1,cosph2,sinph2 !maxsingle
+      real(kind=8),dimension(ndouble_nucl+1,ndouble_nucl+1) :: cosph1ph2,sinph1ph2 !maxdouble,maxdouble
+      logical :: lprn=.true., lprn1=.false.
+!el local variables
+      integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m
+      real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai
+      real(kind=8) :: aux,etheta_nucl,ccl,ssl,scl,csl,ethetacnstr
+! local variables for constrains
+      real(kind=8) :: difi,thetiii
+       integer itheta
+      etheta_nucl=0.0D0
+      print *,"ithet_start",ithet_nucl_start," ithet_end",ithet_nucl_end,nres
+      do i=ithet_nucl_start,ithet_nucl_end
+        if ((itype(i-1,2).eq.ntyp1_molec(2)).or.&
+        (itype(i-2,2).eq.ntyp1_molec(2)).or.     &
+        (itype(i,2).eq.ntyp1_molec(2))) cycle
+        dethetai=0.0d0
+        dephii=0.0d0
+        dephii1=0.0d0
+        theti2=0.5d0*theta(i)
+        ityp2=ithetyp_nucl(itype(i-1,2))
+        do k=1,nntheterm_nucl
+          coskt(k)=dcos(k*theti2)
+          sinkt(k)=dsin(k*theti2)
+        enddo
+        if (i.gt.3 .and. itype(i-2,2).ne.ntyp1_molec(2)) then
+#ifdef OSF
+          phii=phi(i)
+          if (phii.ne.phii) phii=150.0
+#else
+          phii=phi(i)
+#endif
+          ityp1=ithetyp_nucl(itype(i-2,2))
+          do k=1,nsingle_nucl
+            cosph1(k)=dcos(k*phii)
+            sinph1(k)=dsin(k*phii)
+          enddo
+        else
+          phii=0.0d0
+          ityp1=nthetyp_nucl+1
+          do k=1,nsingle_nucl
+            cosph1(k)=0.0d0
+            sinph1(k)=0.0d0
+          enddo
+        endif
+
+        if (i.lt.nres .and. itype(i,2).ne.ntyp1_molec(2)) then
+#ifdef OSF
+          phii1=phi(i+1)
+          if (phii1.ne.phii1) phii1=150.0
+          phii1=pinorm(phii1)
+#else
+          phii1=phi(i+1)
+#endif
+          ityp3=ithetyp_nucl(itype(i,2))
+          do k=1,nsingle_nucl
+            cosph2(k)=dcos(k*phii1)
+            sinph2(k)=dsin(k*phii1)
+          enddo
+        else
+          phii1=0.0d0
+          ityp3=nthetyp_nucl+1
+          do k=1,nsingle_nucl
+            cosph2(k)=0.0d0
+            sinph2(k)=0.0d0
+          enddo
+        endif
+        ethetai=aa0thet_nucl(ityp1,ityp2,ityp3)
+        do k=1,ndouble_nucl
+          do l=1,k-1
+            ccl=cosph1(l)*cosph2(k-l)
+            ssl=sinph1(l)*sinph2(k-l)
+            scl=sinph1(l)*cosph2(k-l)
+            csl=cosph1(l)*sinph2(k-l)
+            cosph1ph2(l,k)=ccl-ssl
+            cosph1ph2(k,l)=ccl+ssl
+            sinph1ph2(l,k)=scl+csl
+            sinph1ph2(k,l)=scl-csl
+          enddo
+        enddo
+        if (lprn) then
+        write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,&
+         " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1
+        write (iout,*) "coskt and sinkt",nntheterm_nucl
+        do k=1,nntheterm_nucl
+          write (iout,*) k,coskt(k),sinkt(k)
+        enddo
+        endif
+        do k=1,ntheterm_nucl
+          ethetai=ethetai+aathet_nucl(k,ityp1,ityp2,ityp3)*sinkt(k)
+          dethetai=dethetai+0.5d0*k*aathet_nucl(k,ityp1,ityp2,ityp3)&
+           *coskt(k)
+          if (lprn)&
+         write (iout,*) "k",k," aathet",aathet_nucl(k,ityp1,ityp2,ityp3),&
+          " ethetai",ethetai
+        enddo
+        if (lprn) then
+        write (iout,*) "cosph and sinph"
+        do k=1,nsingle_nucl
+          write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k)
+        enddo
+        write (iout,*) "cosph1ph2 and sinph2ph2"
+        do k=2,ndouble_nucl
+          do l=1,k-1
+            write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l),&
+              sinph1ph2(l,k),sinph1ph2(k,l)
+          enddo
+        enddo
+        write(iout,*) "ethetai",ethetai
+        endif
+        do m=1,ntheterm2_nucl
+          do k=1,nsingle_nucl
+            aux=bbthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)&
+              +ccthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k)&
+              +ddthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)&
+              +eethet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k)
+            ethetai=ethetai+sinkt(m)*aux
+            dethetai=dethetai+0.5d0*m*aux*coskt(m)
+            dephii=dephii+k*sinkt(m)*(&
+               ccthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)-&
+               bbthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k))
+            dephii1=dephii1+k*sinkt(m)*(&
+               eethet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)-&
+               ddthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k))
+            if (lprn) &
+           write (iout,*) "m",m," k",k," bbthet",&
+              bbthet_nucl(k,m,ityp1,ityp2,ityp3)," ccthet",&
+              ccthet_nucl(k,m,ityp1,ityp2,ityp3)," ddthet",&
+              ddthet_nucl(k,m,ityp1,ityp2,ityp3)," eethet",&
+              eethet_nucl(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+          enddo
+        enddo
+        if (lprn) &
+        write(iout,*) "ethetai",ethetai
+        do m=1,ntheterm3_nucl
+          do k=2,ndouble_nucl
+            do l=1,k-1
+              aux=ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+&
+                 ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+&
+                 ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+&
+                 ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
+              ethetai=ethetai+sinkt(m)*aux
+              dethetai=dethetai+0.5d0*m*coskt(m)*aux
+              dephii=dephii+l*sinkt(m)*(&
+                -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-&
+                 ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+&
+                 ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+&
+                 ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+              dephii1=dephii1+(k-l)*sinkt(m)*( &
+                -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+&
+                 ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+&
+                 ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-&
+                 ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+              if (lprn) then
+              write (iout,*) "m",m," k",k," l",l," ffthet", &
+                 ffthet_nucl(l,k,m,ityp1,ityp2,ityp3), &
+                 ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ggthet",&
+                 ggthet_nucl(l,k,m,ityp1,ityp2,ityp3),&
+                 ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+              write (iout,*) cosph1ph2(l,k)*sinkt(m), &
+                 cosph1ph2(k,l)*sinkt(m),&
+                 sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
+              endif
+            enddo
+          enddo
+        enddo
+10      continue
+        if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') &
+        i,theta(i)*rad2deg,phii*rad2deg, &
+        phii1*rad2deg,ethetai
+        etheta_nucl=etheta_nucl+ethetai
+        print *,i,"partial sum",etheta_nucl
+        if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang_nucl*dephii
+        if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang_nucl*dephii1
+        gloc(nphi+i-2,icg)=wang_nucl*dethetai
+      enddo
+      return
+      end subroutine ebend_nucl
+
  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
        end module energy
  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
        end module energy