From: Adam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Date: Wed, 17 Dec 2014 09:48:55 +0000 (+0100)
Subject: Merge branch 'master' of mmka:unres
X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=78c6b598700d2c701af6532afb20ebc905a0b8ef;hp=-c;p=unres.git

Merge branch 'master' of mmka:unres

Conflicts:
	.gitignore
	bin/cluster/unres_clustMD-mult_MPICH-GAB.exe
	bin/unres/MD/unres_gfortran_single_GAB.exe
	bin/unres/MD/unres_ifort_single_GAB.exe
	bin/unres/MINIM/unres_ifort_MIN_single_E0LL2Y.exe
	bin/unres/MINIM/unres_ifort_MIN_single_GAB.exe
	bin/unres_clustMD_MPI-oldparm
	bin/wham/wham_multparm-ham_rep-oldparm
	source/cluster/wham/src/COMMON.SCCOR
	source/cluster/wham/src/Makefile
	source/unres/src_MD-M/Makefile
	source/unres/src_MD-M/energy_p_new_barrier.F
	source/unres/src_MD-M/stochfric.F
	source/unres/src_MD/cinfo.f
	source/wham/src-M/Makefile
	source/wham/src/DIMENSIONS.FREE
---

78c6b598700d2c701af6532afb20ebc905a0b8ef
diff --combined .gitignore
index 442c6c4,2d9541f..0604082
--- a/.gitignore
+++ b/.gitignore
@@@ -12,8 -12,7 +12,9 @@@ cinfo.
  # ignore build dir
  build/
  build2/
 +build_*/
 +period_*/
+ 
  # latex files in documentation 
  doc/*/*.aux
  doc/*/*.log
@@@ -24,7 -23,6 +25,7 @@@ gradcheck
  mapcheck/
  run/
  sympcheck/
 -compinfo
 -DIL/
  bin/unres/MD/unres_ifort_MPICH_GAB_czyt.exe
 +bin/unres/MD-M/unres_Tc_procor_newparm_em64-D-symetr.exe
 +DIL/
 +compinfo
diff --combined CMakeLists.txt
index 9b6b6b8,5b824dd..89bd8ef
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@@ -6,8 -6,8 +6,8 @@@ cmake_minimum_required(VERSION 2.8
  project(UNRESPACK Fortran C)
  
  set(UNRES_MAJOR 3)
- set(UNRES_MINOR 1)
- set(UNRES_PATCH 0)
+ set(UNRES_MINOR 2)
+ set(UNRES_PATCH 1)
  set(UNRES_VERSION ${UNRES_MAJOR}.${UNRES_MINOR}.${UNRES_PATCH})
  
  #======================================
@@@ -63,6 -63,13 +63,13 @@@ MACRO (CINFO_FORMAT FN VN VD
  	file(APPEND ${FN} "       write(iout,*)'${STR}'\n")
      endif(SUMA GREATER 50)
  ENDMACRO (CINFO_FORMAT)
+ 
+ # Some MPI wrappers pass double include paths
+ # This macro fixes broken by semicolon occurence in path
+ MACRO (FIX_DBL_INCLUDE RESULT)
+   string(REPLACE ";" " -I" ${RESULT} "${${RESULT}}")
+ ENDMACRO (FIX_DBL_INCLUDE)
+ 
  #======================================
  # CTest stuff
  #======================================A
@@@ -70,8 -77,6 +77,6 @@@
  include(CTest)
  enable_testing()
   
- # Set makefile verbose on
- set( CMAKE_VERBOSE_MAKEFILE 1 )
  
  #======================================
  # Fortran compilers stuff
@@@ -88,10 -93,23 +93,23 @@@ SET(CMAKE_Fortran_COMPILE_OBJECT "<CMAK
     
  # make sure that the default is a RELEASE
  if (NOT CMAKE_BUILD_TYPE)
-   set (CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: None Debug Release." FORCE)
-   set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "None" "Debug" "Release" )
+   set (CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: None Debug Release RelWithDebInfo." FORCE)
+   set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "None" "Debug" "Release" "RelWithDebInfo")
  endif (NOT CMAKE_BUILD_TYPE)
-    
+ 
+ 
+ if (CMAKE_BUILD_TYPE STREQUAL "Release")
+   # Set makefile verbosity off for Release builds
+   set( CMAKE_VERBOSE_MAKEFILE 0 )
+ else()
+   # Set makefile verbosity on for other builds
+   set( CMAKE_VERBOSE_MAKEFILE 1 )
+ endif (CMAKE_BUILD_TYPE STREQUAL "Release")
+ 
+ # Default Install Path
+ 
+ set(CMAKE_INSTALL_PREFIX "${CMAKE_SOURCE_DIR}/bin" CACHE PATH "Binary install directory " FORCE)
+ 
  #=======================================  
  # Set the varous build variables 
  #=======================================
@@@ -120,54 -138,16 +138,16 @@@ option(UNRES_NA_MMCE "Kompilujemy na mm
  # MPI stuff
  #=================================
  
- # Note for the future - use finde package to get MPI 
- find_package(MPI)
- 
- #if(MPI_LIBRARY)
- 
- #MPI_INCLUDE_PATH
- 
- 
- if(MPIF_LOCAL_DIR)
-   find_library(MPIF_LIBRARY NAMES libmpich.a  NO_DEFAULT_PATH  PATHS  ${MPIF_LOCAL_DIR}/lib)
-   find_path( MPIF_INCLUDE_DIRECTORIES  NAMES mpif.h  NO_DEFAULT_PATH  PATHS  ${MPIF_LOCAL_DIR}/include  )
- else(MPIF_LOCAL_DIR)
-   find_library(MPIF_LIBRARY NAMES mpi mpich PATHS 
- 	${MPI_LIBRARY}
- 	${MPI_LIBRARY}/../
- 	${MPI_EXTRA_LIBRARY}
- 	/users/local/mpi64/mpich-1.2.7p1/lib 
- 	/usr/lib
-         /usr/local/lib
-         /usr/local/mpi/lib
-   )
-   find_path( MPIF_INCLUDE_DIRECTORIES NAMES  mpif.h PATHS
- 	${MPI_INCLUDE_PATH}
- 	/users/local/mpi64/mpich-1.2.7p1/include
- 	/usr/include
- 	/usr/local/include   
- 	/usr/include/mpi
- 	/usr/local/mpi/include
-   )
- endif(MPIF_LOCAL_DIR)
- 
- set( MPIF_LIBRARIES  ${MPIF_LIBRARY})
- 
- if ( MPIF_INCLUDE_DIRECTORIES )
-   set( MPIF_FOUND TRUE )
+ # Note for the future - use find package to get MPI 
+ find_package(MPI QUIET)
+ 
+ if (MPI_Fortran_FOUND)
    message("MPI found")
- else ( MPIF_INCLUDE_DIRECTORIES )
-   set( MPIF_FOUND FALSE )
+   FIX_DBL_INCLUDE(MPI_Fortran_INCLUDE_PATH)
+ else()
    message("MPI not found - disabling MPI compile flags ")
    set ( UNRES_WITH_MPI "OFF")
- endif ( MPIF_INCLUDE_DIRECTORIES )
- 
- if (MPIF_FOUND)
-   message("MPIF_LIBRARIES=${MPIF_LIBRARY}")
-   message("MPIF_INCLUDE_DIRECTORIES=${MPIF_INCLUDE_DIRECTORIES}" )
- endif(MPIF_FOUND) 
- 
- 
+ endif(MPI_Fortran_FOUND)	
  
  #======================================
  # Detect system architecture
@@@ -188,29 -168,22 +168,21 @@@ message("Detected ${architektura}-bit a
  # used by unres/src_MIN
  find_package (Threads)
  
- 
- 
- #=======================================
- #  Create diractories for build targets
- #=======================================
- #execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_BINARY_DIR}/bin/unres/MD )
- 
  #=======================================
  # Add source files
  #=======================================     
  
 -
  add_subdirectory(source/lib/xdrf)
  
 -
  if(UNRES_NA_MMCE)
  
    if(UNRES_WITH_MPI)
      # Brak MPI dla gfortrana, wiec tylko na ifort sie skompiluje
      if (Fortran_COMPILER_NAME STREQUAL "ifort")
        add_subdirectory(source/unres/src_MD)
 +      add_subdirectory(source/unres/src_MD_DFA)
        add_subdirectory(source/unres/src_MD-M)
        add_subdirectory(source/unres/src_CSA)
-       add_subdirectory(source/unres/src_CSA_DiL)
        add_subdirectory(source/cluster/wham/src)
        add_subdirectory(source/cluster/wham/src-M)
      endif (Fortran_COMPILER_NAME STREQUAL "ifort")
@@@ -229,9 -202,7 +201,8 @@@ else(
    add_subdirectory(source/unres/src_MD)
    if(UNRES_WITH_MPI)
      add_subdirectory(source/unres/src_MD-M)
 +    add_subdirectory(source/unres/src_MD_DFA)
      add_subdirectory(source/unres/src_CSA)
-     add_subdirectory(source/unres/src_CSA_DiL)
      add_subdirectory(source/wham/src)
      add_subdirectory(source/wham/src-M)
      add_subdirectory(source/cluster/wham/src)
diff --combined bin/unres/MINIM/unres_ifort_MIN_single_GAB.exe
index 6e1ebca,8fbedcf..0000000
deleted file mode 100755,100755
Binary files differ
diff --combined source/unres/src_MD-M/CMakeLists.txt
index 3884bc1,ecd2b5a..de40196
--- a/source/unres/src_MD-M/CMakeLists.txt
+++ b/source/unres/src_MD-M/CMakeLists.txt
@@@ -5,11 -5,6 +5,6 @@@
  enable_language (Fortran)
  
  #================================
- # build the xdrf library 
- #================================ 
- #add_subdirectory(xdrf)
- 
- #================================
  # Set source file lists
  #================================
  set(UNRES_MDM_SRC0 
@@@ -71,6 -66,7 +66,7 @@@
  	permut.F
  	pinorm.f 
  	printmat.f 
+ 	prng_32.F
  	q_measure.F 
  	ran.f
  	randgens.f 
@@@ -93,17 -89,9 +89,9 @@@
  	timing.F
  	together.F
  	unres.F
+ 	ssMD.F
  )
  
- if (Fortran_COMPILER_NAME STREQUAL "ifort")
-   set(UNRES_MDM_SRC0 ${UNRES_MDM_SRC0} prng.f ) 
- elseif(Fortran_COMPILER_NAME STREQUAL "mpif90")
-   set(UNRES_MDM_SRC0 ${UNRES_MDM_SRC0} prng.f )
- else()
-   set(UNRES_MDM_SRC0 ${UNRES_MDM_SRC0} prng_32.F )
- endif (Fortran_COMPILER_NAME STREQUAL "ifort")
- 
- 
  set(UNRES_MDM_SRC3 energy_p_new_barrier.F energy_p_new-sep_barrier.F gradient_p.F )
  
  set(UNRES_MDM_PP_SRC
@@@ -141,6 -129,7 +129,7 @@@
  	newconf.f 
  	parmread.F 
  	permut.F
+ 	prng_32.F
  	q_measure1.F 
  	q_measure3.F 
  	q_measure.F
@@@ -187,10 -176,10 +176,10 @@@ endif (Fortran_COMPILER_NAME STREQUAL "
  
  # Add MPI compiler flags
  if(UNRES_WITH_MPI)
-   set(FFLAGS0 "${FFLAGS0} -I${MPIF_INCLUDE_DIRECTORIES}")
-   set(FFLAGS1 "${FFLAGS1} -I${MPIF_INCLUDE_DIRECTORIES}")
-   set(FFLAGS2 "${FFLAGS2} -I${MPIF_INCLUDE_DIRECTORIES}")
-   set(FFLAGS3 "${FFLAGS3} -I${MPIF_INCLUDE_DIRECTORIES}")
+   set(FFLAGS0 "${FFLAGS0} -I${MPI_Fortran_INCLUDE_PATH}")
+   set(FFLAGS1 "${FFLAGS1} -I${MPI_Fortran_INCLUDE_PATH}")
+   set(FFLAGS2 "${FFLAGS2} -I${MPI_Fortran_INCLUDE_PATH}")
+   set(FFLAGS3 "${FFLAGS3} -I${MPI_Fortran_INCLUDE_PATH}")
  endif(UNRES_WITH_MPI)
  
  set_property(SOURCE ${UNRES_MDM_SRC0} APPEND PROPERTY COMPILE_FLAGS ${FFLAGS0} )
@@@ -203,14 -192,14 +192,14 @@@ set_property(SOURCE ${UNRES_MDM_SRC3} P
  #=========================================
  if(UNRES_MD_FF STREQUAL "GAB" )
    # set preprocesor flags   
 -  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC" )
 +  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC  -DSCCORPDB -DNEWCORR" )
  
  #=========================================
  #  Settings for E0LL2Y force field
  #=========================================
  elseif(UNRES_MD_FF STREQUAL "E0LL2Y")
    # set preprocesor flags   
 -  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0" )
 +  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DSCCORPDB -DNEWCORR" )
  endif(UNRES_MD_FF STREQUAL "GAB")
  
  
@@@ -235,9 -224,6 +224,9 @@@ elseif (Fortran_COMPILER_NAME STREQUAL 
  elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
    # Add old gfortran flags
    set(CPPFLAGS "${CPPFLAGS} -DG77") 
 +else (Fortran_COMPILER_NAME STREQUAL "ifort")
 +  # Default preprocessor flags
 +  set(CPPFLAGS "${CPPFLAGS} -DPGI")
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
  
  
@@@ -249,7 -235,9 +238,9 @@@ if (UNRES_WITH_MPI
  endif(UNRES_WITH_MPI)
  
  
+ #=========================================
  # add 64-bit specific preprocessor flags
+ #=========================================
  if (architektura STREQUAL "64")
    set(CPPFLAGS "${CPPFLAGS} -DAMD64")
  endif (architektura STREQUAL "64")
@@@ -268,10 -256,10 +259,10 @@@ set_property(SOURCE proc_proc.c PROPERT
  #========================================
  if(UNRES_WITH_MPI) 
    # binary with mpi
-   set(UNRES_BIN "unres_${Fortran_COMPILER_NAME}_MPICH_${UNRES_FF}.exe")
+   set(UNRES_BIN "unresMD-M_${Fortran_COMPILER_NAME}_MPICH_${UNRES_MD_FF}.exe")
  else(UNRES_WITH_MPI)
    # binary without mpi
-   set(UNRES_BIN "unres_${Fortran_COMPILER_NAME}_single_${UNRES_FF}.exe")
+   set(UNRES_BIN "unresMD-M_${Fortran_COMPILER_NAME}_single_${UNRES_MD_FF}.exe")
  endif(UNRES_WITH_MPI)  
  
  #=========================================
@@@ -316,8 -304,7 +307,7 @@@ set(UNRES_MDM_SRCS ${UNRES_MDM_SRC0} ${
  #=========================================
  add_executable(UNRES_BIN-MD-M ${UNRES_MDM_SRCS} )
  set_target_properties(UNRES_BIN-MD-M PROPERTIES OUTPUT_NAME ${UNRES_BIN})
- 
- #set_property(TARGET ${UNRES_BIN} PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin/unres/MD )
+ set_property(TARGET UNRES_BIN-MD-M PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  #add_dependencies (${UNRES_BIN} ${UNRES_XDRFLIB})
  
  #=========================================
@@@ -325,11 -312,16 +315,16 @@@
  #=========================================
  # link MPI library (libmpich.a)  
  if(UNRES_WITH_MPI)
-   target_link_libraries( UNRES_BIN-MD-M ${MPIF_LIBRARIES} )
+   target_link_libraries( UNRES_BIN-MD-M ${MPI_Fortran_LIBRARIES} )
  endif(UNRES_WITH_MPI)
  # link libxdrf.a 
  #message("UNRES_XDRFLIB=${UNRES_XDRFLIB}")
  target_link_libraries( UNRES_BIN-MD-M xdrf )
+ 
+ #=========================================
+ # Install Path
+ #=========================================
+ install(TARGETS UNRES_BIN-MD-M DESTINATION ${CMAKE_INSTALL_PREFIX})
  
  #=========================================
  # TESTS 
diff --combined source/unres/src_MD-M/MD_A-MTS.F
index 1144934,56b3ea8..1191c18
--- a/source/unres/src_MD-M/MD_A-MTS.F
+++ b/source/unres/src_MD-M/MD_A-MTS.F
@@@ -209,7 -209,7 +209,7 @@@ c Variable time step algorithm
            enddo
          enddo
          do i=nnt,nct
 -          if (itype(i).ne.10 .and. itype(i).ne.21) then
 +          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
              do j=1,3
                ind=ind+1
                v_work(ind)=d_t(j,i+nres)
@@@ -955,7 -955,7 +955,7 @@@ c  forces)
          enddo
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            inres=i+nres
            do j=1,3
              d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
@@@ -1005,7 -1005,7 +1005,7 @@@ c Applying velocity Verlet algorithm - 
          enddo
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            inres=i+nres
            do j=1,3    
              adt=d_a_old(j,inres)*d_time
@@@ -1049,7 -1049,7 +1049,7 @@@ c  Step 2 of the velocity Verlet algori
          enddo
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            inres=i+nres
            do j=1,3
              d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
@@@ -1087,12 -1087,25 +1087,25 @@@ c Applying velocity Verlet algorithm - 
  c
  c Compute friction and stochastic forces
  c
+ #ifdef MPI
        time00=MPI_Wtime()
+ #else
+       time00=tcpu()
+ #endif
        call friction_force
+ #ifdef MPI
        time_fric=time_fric+MPI_Wtime()-time00
        time00=MPI_Wtime()
+ #else
+       time_fric=time_fric+tcpu()-time00
+       time00=tcpu()
+ #endif
        call stochastic_force(stochforcvec) 
+ #ifdef MPI
        time_stoch=time_stoch+MPI_Wtime()-time00
+ #else
+       time_stoch=time_stoch+tcpu()-time00
+ #endif
  c
  c Compute the acceleration due to friction forces (d_af_work) and stochastic
  c forces (d_as_work)
@@@ -1148,7 -1161,7 +1161,7 @@@
          ind=ind+3
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            inres=i+nres
            do j=1,3    
              adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
@@@ -1213,7 -1226,7 +1226,7 @@@
          ind=ind+3
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            inres=i+nres
            do j=1,3
              d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
@@@ -1260,7 -1273,6 +1273,7 @@@ c            accel(j)=aux(j)+0.5d0*(d_a
  c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
              if (dabs(accel(j)).gt.dabs(accel_old(j))) then
                dacc=dabs(accel(j)-accel_old(j))
 +c              write (iout,*) i,dacc
                if (dacc.gt.amax) amax=dacc
              endif
            enddo
@@@ -1279,7 -1291,7 +1292,7 @@@ c        accel(j)=aux(j
          enddo
        endif
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            do j=1,3 
  c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
              accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
@@@ -1290,7 -1302,6 +1303,7 @@@
  c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
            if (dabs(accel(j)).gt.dabs(accel_old(j))) then
              dacc=dabs(accel(j)-accel_old(j))
 +c            write (iout,*) "side-chain",i,dacc
              if (dacc.gt.amax) amax=dacc
            endif
          enddo
@@@ -1333,7 -1344,7 +1346,7 @@@ c            write (iout,*) "back",i,j,
            enddo
          endif
  c Side chains
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            do j=1,3 
              epdriftij=
       &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
@@@ -1380,7 -1391,7 +1393,7 @@@ c      write(iout,*) "fact", fac
          enddo
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            inres=i+nres
            do j=1,3
              d_t(j,inres)=fact*d_t(j,inres)
@@@ -1435,8 -1446,7 +1448,8 @@@ c if the friction coefficients do not d
            stdforcp(i)=stdfp*dsqrt(gamp)
          enddo
          do i=nnt,nct
 -          stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
 +          stdforcsc(i)=stdfsc(iabs(itype(i)))
 +     &                *dsqrt(gamsc(iabs(itype(i))))
          enddo
        endif
  c Open the pdb file for snapshotshots
@@@ -1510,11 -1520,13 +1523,13 @@@ c        statname=statname(:ilen(statna
         if (restart1file) then
           if (me.eq.king)
       &     inquire(file=mremd_rst_name,exist=file_exist)
+ #ifdef MPI
             write (*,*) me," Before broadcast: file_exist",file_exist
           call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
       &          IERR)
           write (*,*) me," After broadcast: file_exist",file_exist
  c        inquire(file=mremd_rst_name,exist=file_exist)
+ #endif
          if(me.eq.king.or..not.out1file)
       &   write(iout,*) "Initial state read by master and distributed"
         else
@@@ -1827,7 -1839,7 +1842,7 @@@ c Transfer to the d_t vecto
          enddo
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            do j=1,3
              ind=ind+1
              d_t(j,i+nres)=d_t_work(ind)
@@@ -2056,7 -2068,7 +2071,7 @@@ c      endd
          ind=ind+3
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            do j=1,3
              dc_work(ind+j)=dc_old(j,i+nres)
              d_t_work(ind+j)=d_t_old(j,i+nres)
@@@ -2165,7 -2177,7 +2180,7 @@@
          ind=ind+3
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            inres=i+nres
            do j=1,3
              d_t(j,inres)=d_t_work(ind+j)
@@@ -2322,7 -2334,7 +2337,7 @@@ c      endd
          ind=ind+3
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            do j=1,3
              dc_work(ind+j)=dc_old(j,i+nres)
              d_t_work(ind+j)=d_t_old(j,i+nres)
@@@ -2371,7 -2383,7 +2386,7 @@@
          ind=ind+3
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            inres=i+nres
            do j=1,3
              dc(j,inres)=dc_work(ind+j)
@@@ -2432,7 -2444,7 +2447,7 @@@ c          ddt2=ddt2+vrand_mat2(i,j)*st
          ind=ind+3
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            inres=i+nres
            do j=1,3
              d_t(j,inres)=d_t_work(ind+j)
diff --combined source/unres/src_MD-M/energy_p_new_barrier.F
index 27c5f1d,2097265..6976071
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@@ -121,6 -121,11 +121,11 @@@
  C Calculate electrostatic (H-bonding) energy of the main chain.
  C
    107 continue
+ cmc
+ cmc Sep-06: egb takes care of dynamic ss bonds too
+ cmc
+ c      if (dyn_ss) call dyn_set_nss
+ 
  c      print *,"Processor",myrank," computed USCSC"
  #ifdef TIMING
        time01=MPI_Wtime() 
@@@ -296,10 -301,9 +301,11 @@@
        energia(17)=estr
        energia(20)=Uconst+Uconst_back
        energia(21)=esccor
 +c    Here are the energies showed per procesor if the are more processors 
 +c    per molecule then we sum it up in sum_energy subroutine 
  c      print *," Processor",myrank," calls SUM_ENERGY"
        call sum_energy(energia,.true.)
+       if (dyn_ss) call dyn_set_nss
  c      print *," Processor",myrank," left SUM_ENERGY"
  #ifdef TIMING
        time_sumene=time_sumene+MPI_Wtime()-time00
@@@ -435,9 -439,9 +441,9 @@@ cMS$ATTRIBUTES C ::  proc_pro
  #endif
  #ifdef MPI
        include 'mpif.h'
+ #endif
        double precision gradbufc(3,maxres),gradbufx(3,maxres),
 -     &  glocbuf(4*maxres),gradbufc_sum(3,maxres)
 +     &  glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
- #endif
        include 'COMMON.SETUP'
        include 'COMMON.IOUNITS'
        include 'COMMON.FFIELD'
@@@ -449,7 -453,6 +455,7 @@@
        include 'COMMON.CONTROL'
        include 'COMMON.TIME1'
        include 'COMMON.MAXGRAD'
 +      include 'COMMON.SCCOR'
  #ifdef TIMING
        time01=MPI_Wtime()
  #endif
@@@ -692,6 -695,7 +698,6 @@@ c      endd
       &   +wturn3*gel_loc_turn3(i)
       &   +wturn6*gel_loc_turn6(i)
       &   +wel_loc*gel_loc_loc(i)
 -     &   +wsccor*gsccor_loc(i)
        enddo
  #ifdef DEBUG
        write (iout,*) "gloc after adding corr"
@@@ -710,21 -714,6 +716,21 @@@
          do i=1,4*nres
            glocbuf(i)=gloc(i,icg)
          enddo
 +#define DEBUG
 +#ifdef DEBUG
 +      write (iout,*) "gloc_sc before reduce"
 +      do i=1,nres
 +       do j=1,1
 +        write (iout,*) i,j,gloc_sc(j,i,icg)
 +       enddo
 +      enddo
 +#endif
 +#undef DEBUG
 +        do i=1,nres
 +         do j=1,3
 +          gloc_scbuf(j,i)=gloc_sc(j,i,icg)
 +         enddo
 +        enddo
          time00=MPI_Wtime()
          call MPI_Barrier(FG_COMM,IERR)
          time_barrier_g=time_barrier_g+MPI_Wtime()-time00
@@@ -736,19 -725,6 +742,19 @@@
          call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,
       &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
          time_reduce=time_reduce+MPI_Wtime()-time00
 +        call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,
 +     &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
 +        time_reduce=time_reduce+MPI_Wtime()-time00
 +#define DEBUG
 +#ifdef DEBUG
 +      write (iout,*) "gloc_sc after reduce"
 +      do i=1,nres
 +       do j=1,1
 +        write (iout,*) i,j,gloc_sc(j,i,icg)
 +       enddo
 +      enddo
 +#endif
 +#undef DEBUG
  #ifdef DEBUG
        write (iout,*) "gloc after reduce"
        do i=1,4*nres
@@@ -1055,9 -1031,9 +1061,9 @@@
  c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        if (itypi.eq.21) cycle
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        if (itypi.eq.ntyp1) cycle
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -1070,8 -1046,8 +1076,8 @@@
  cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  cd   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
 -            itypj=itype(j)
 -            if (itypj.eq.21) cycle
 +            itypj=iabs(itype(j)) 
 +            if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              zj=c(3,nres+j)-zi
@@@ -1208,9 -1184,9 +1214,9 @@@
  c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        if (itypi.eq.21) cycle
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        if (itypi.eq.ntyp1) cycle
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -1219,8 -1195,8 +1225,8 @@@ C Calculate SC interaction energy
  C
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
 -            itypj=itype(j)
 -            if (itypj.eq.21) cycle
 +            itypj=iabs(itype(j))
 +            if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              zj=c(3,nres+j)-zi
@@@ -1301,9 -1277,9 +1307,9 @@@ c     els
  c     endif
        ind=0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        if (itypi.eq.21) cycle
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        if (itypi.eq.ntyp1) cycle
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -1318,8 -1294,8 +1324,8 @@@
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
              ind=ind+1
 -            itypj=itype(j)
 -            if (itypj.eq.21) cycle
 +            itypj=iabs(itype(j))
 +            if (itypj.eq.ntyp1) cycle
  c            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              chi1=chi(itypi,itypj)
@@@ -1412,6 -1388,7 +1418,7 @@@
        include 'COMMON.IOUNITS'
        include 'COMMON.CALC'
        include 'COMMON.CONTROL'
+       include 'COMMON.SBRIDGE'
        logical lprn
        evdw=0.0D0
  ccccc      energy_dec=.false.
@@@ -1421,9 -1398,9 +1428,9 @@@ c     print *,'Entering EGB nnt=',nnt,
  c     if (icall.eq.0) lprn=.false.
        ind=0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        if (itypi.eq.21) cycle
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        if (itypi.eq.ntyp1) cycle
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -1439,9 -1416,15 +1446,15 @@@ C Calculate SC interaction energy
  C
          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'
+             ELSE
              ind=ind+1
 -            itypj=itype(j)
 -            if (itypj.eq.21) cycle
 +            itypj=iabs(itype(j))
 +            if (itypj.eq.ntyp1) cycle
  c            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
  c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
@@@ -1533,6 -1516,7 +1546,7 @@@ C Calculate the radial part of the grad
              gg(3)=zj*fac
  C Calculate angular part of the gradient.
              call sc_grad
+             ENDIF    ! dyn_ss            
            enddo      ! j
          enddo        ! iint
        enddo          ! i
@@@ -1566,9 -1550,9 +1580,9 @@@ c     print *,'Entering EGB nnt=',nnt,
  c     if (icall.eq.0) lprn=.true.
        ind=0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        if (itypi.eq.21) cycle
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        if (itypi.eq.ntyp1) cycle
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -1583,8 -1567,8 +1597,8 @@@
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
              ind=ind+1
 -            itypj=itype(j)
 -            if (itypj.eq.21) cycle
 +            itypj=iabs(itype(j))
 +            if (itypj.eq.ntyp1) cycle
  c            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              sig0ij=sigma(itypi,itypj)
@@@ -1814,9 -1798,9 +1828,9 @@@
  cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
        evdw=0.0D0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        if (itypi.eq.21) cycle
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        if (itypi.eq.ntyp1) cycle
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -1827,8 -1811,8 +1841,8 @@@
  cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  cd   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
 -            itypj=itype(j)
 -            if (itypj.eq.21) cycle
 +            itypj=iabs(itype(j))
 +            if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              zj=c(3,nres+j)-zi
@@@ -1896,7 -1880,7 +1910,7 @@@ cd      write(iout,*) 'In EELEC_soft_sp
        eello_turn4=0.0d0
        ind=0
        do i=iatel_s,iatel_e
 -        if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
 +        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@@ -1906,7 -1890,7 +1920,7 @@@
          num_conti=0
  c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
          do j=ielstart(i),ielend(i)
 -          if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle
 +          if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
            ind=ind+1
            iteli=itel(i)
            itelj=itel(j)
@@@ -2252,78 -2236,11 +2266,78 @@@
  C Compute the virtual-bond-torsional-angle dependent quantities needed
  C to calculate the el-loc multibody terms of various order.
  C
 +c      write(iout,*) 'nphi=',nphi,nres
 +#ifdef PARMAT
 +      do i=ivec_start+2,ivec_end+2
 +#else
 +      do i=3,nres+1
 +#endif
 +#ifdef NEWCORR
 +        if (i.gt. nnt+2 .and. i.lt.nct+2) then
 +          iti = itortyp(itype(i-2))
 +        else
 +          iti=ntortyp+1
 +        endif
 +c        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))
 +        else
 +          iti1=ntortyp+1
 +        endif
 +c        write(iout,*),i
 +        b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0)
 +     &           +bnew1(2,1,iti)*dsin(theta(i-1))
 +     &           +bnew1(3,1,iti)*dcos(theta(i-1)/2.0)
 +        gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
 +     &             +bnew1(2,1,iti)*dcos(theta(i-1))
 +     &             -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
 +c     &           +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i))
 +c     &*(cos(theta(i)/2.0)
 +        b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0)
 +     &           +bnew2(2,1,iti)*dsin(theta(i-1))
 +     &           +bnew2(3,1,iti)*dcos(theta(i-1)/2.0)
 +c     &           +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i))
 +c     &*(cos(theta(i)/2.0)
 +        gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
 +     &             +bnew2(2,1,iti)*dcos(theta(i-1))
 +     &             -bnew2(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
 +c        if (ggb1(1,i).eq.0.0d0) then
 +c        write(iout,*) 'i=',i,ggb1(1,i),
 +c     &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0,
 +c     &bnew1(2,1,iti)*cos(theta(i)),
 +c     &bnew1(3,1,iti)*sin(theta(i)/2.0)/2.0
 +c        endif
 +        b1(2,i-2)=bnew1(1,2,iti)
 +        gtb1(2,i-2)=0.0
 +        b2(2,i-2)=bnew2(1,2,iti)
 +        gtb2(2,i-2)=0.0
 +        EE(1,1,i-2)=eenew(1,iti)*dcos(theta(i-1))
 +        EE(1,2,i-2)=eeold(1,2,iti)
 +        EE(2,1,i-2)=eeold(2,1,iti)
 +        EE(2,2,i-2)=eeold(2,2,iti)
 +        gtEE(1,1,i-2)=-eenew(1,iti)*dsin(theta(i-1))
 +        gtEE(1,2,i-2)=0.0d0
 +        gtEE(2,2,i-2)=0.0d0
 +        gtEE(2,1,i-2)=0.0d0
 +c        EE(2,2,iti)=0.0d0
 +c        EE(1,2,iti)=0.5d0*eenew(1,iti)
 +c        EE(2,1,iti)=0.5d0*eenew(1,iti)
 +c        b1(2,iti)=bnew1(1,2,iti)*sin(alpha(i))*sin(beta(i))
 +c        b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i))
 +       b1tilde(1,i-2)=b1(1,i-2)
 +       b1tilde(2,i-2)=-b1(2,i-2)
 +       b2tilde(1,i-2)=b2(1,i-2)
 +       b2tilde(2,i-2)=-b2(2,i-2)
 +c       write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
 +c       write(iout,*)  'b1=',b1(1,i-2)
 +c       write (iout,*) 'theta=', theta(i-1)
 +       enddo
  #ifdef PARMAT
        do i=ivec_start+2,ivec_end+2
  #else
        do i=3,nres+1
  #endif
 +#endif
          if (i .lt. nres+1) then
            sin1=dsin(phi(i))
            cos1=dcos(phi(i))
@@@ -2408,18 -2325,8 +2422,18 @@@ cd        write (iout,*) 'b2',b2(:,iti
  cd        write (iout,*) 'Ug',Ug(:,:,i-2)
  c        if (i .gt. iatel_s+2) then
          if (i .gt. nnt+2) then
 -          call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2))
 -          call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
 +          call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2))
 +#ifdef NEWCORR
 +          call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2))
 +c          write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj"
 +#endif
 +c          write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti),
 +c     &    EE(1,2,iti),EE(2,2,iti)
 +          call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2))
 +          call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2))
 +c          write(iout,*) "Macierz EUG",
 +c     &    eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2),
 +c     &    eug(2,2,i-2)
            if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) 
       &    then
            call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2))
@@@ -2441,25 -2348,21 +2455,25 @@@
              enddo
            enddo
          endif
 -        call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2))
 -        call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2))
 +        call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2))
 +        call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2))
          do k=1,2
            muder(k,i-2)=Ub2der(k,i-2)
          enddo
  c        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))
 +          if (itype(i-1).le.ntyp) then
 +            iti1 = itortyp(itype(i-1))
 +          else
 +            iti1=ntortyp+1
 +          endif
          else
            iti1=ntortyp+1
          endif
          do k=1,2
 -          mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
 +          mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
          enddo
 -cd        write (iout,*) 'mu ',mu(:,i-2)
 +c        write (iout,*) 'mu ',mu(:,i-2),i-2
  cd        write (iout,*) 'mu1',mu1(:,i-2)
  cd        write (iout,*) 'mu2',mu2(:,i-2)
          if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
@@@ -2470,7 -2373,7 +2484,7 @@@
          call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2))
          call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2))
  C Vectors and matrices dependent on a single virtual-bond dihedral.
 -        call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1))
 +        call matvec2(DD(1,1,iti),b1tilde(1,i-1),auxvec(1))
          call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) 
          call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) 
          call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2))
@@@ -2786,7 -2689,7 +2800,7 @@@
        dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
       &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
        double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
 -     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
 +     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij(4)
        common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
       &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
       &    num_conti,j1,j2
@@@ -2866,8 -2769,8 +2880,8 @@@
  C Loop over i,i+2 and i,i+3 pairs of the peptide groups
  C
        do i=iturn3_start,iturn3_end
 -        if (itype(i).eq.21 .or. itype(i+1).eq.21 
 -     &  .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle
 +        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 +     &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@@ -2883,9 -2786,9 +2897,9 @@@
          num_cont_hb(i)=num_conti
        enddo
        do i=iturn4_start,iturn4_end
 -        if (itype(i).eq.21 .or. itype(i+1).eq.21
 -     &    .or. itype(i+3).eq.21
 -     &    .or. itype(i+4).eq.21) cycle
 +        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 +     &    .or. itype(i+3).eq.ntyp1
 +     &    .or. itype(i+4).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@@ -2896,9 -2799,8 +2910,9 @@@
          ymedi=c(2,i)+0.5d0*dyi
          zmedi=c(3,i)+0.5d0*dzi
          num_conti=num_cont_hb(i)
 +c        write(iout,*) "JESTEM W PETLI"
          call eelecij(i,i+3,ees,evdw1,eel_loc)
 -        if (wturn4.gt.0.0d0 .and. itype(i+2).ne.21) 
 +        if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) 
       &   call eturn4(i,eello_turn4)
          num_cont_hb(i)=num_conti
        enddo   ! i
@@@ -2906,8 -2808,7 +2920,8 @@@
  c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
  c
        do i=iatel_s,iatel_e
 -        if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
 +c       do i=7,7
 +        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@@ -2920,9 -2821,8 +2934,9 @@@
  c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
          num_conti=num_cont_hb(i)
          do j=ielstart(i),ielend(i)
 -c          write (iout,*) i,j,itype(i),itype(j)
 -          if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle
 +c         do j=13,13
 +c          write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
 +          if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
            call eelecij(i,j,ees,evdw1,eel_loc)
          enddo ! j
          num_cont_hb(i)=num_conti
@@@ -2960,8 -2860,7 +2974,8 @@@ C--------------------------------------
        dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
       &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
        double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
 -     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
 +     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
 +     &    gmuij2(4),gmuji2(4)
        common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
       &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
       &    num_conti,j1,j2
@@@ -3025,9 -2924,7 +3039,9 @@@ cd     &      1.0D0/dsqrt(rrmij),evdwij
  cd     &      xmedi,ymedi,zmedi,xj,yj,zj
  
            if (energy_dec) then 
 -              write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
 +              write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') 
 +     &'evdw1',i,j,evdwij
 +     &,iteli,itelj,aaa,evdw1
                write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
            endif
  
@@@ -3178,7 -3075,6 +3192,7 @@@ C   Fourier series in the angles lambda
  C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
  C   are computed for EVERY pair of non-contiguous peptide groups.
  C
 +
            if (j.lt.nres-1) then
              j1=j+1
              j2=j-1
@@@ -3187,20 -3083,10 +3201,20 @@@
              j2=j-2
            endif
            kkk=0
 +          lll=0
            do k=1,2
              do l=1,2
                kkk=kkk+1
                muij(kkk)=mu(k,i)*mu(l,j)
 +c              write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l
 +#ifdef NEWCORR
 +             gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
 +c             write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
 +             gmuij2(kkk)=gUb2(k,i)*mu(l,j)
 +             gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
 +c             write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
 +             gmuji2(kkk)=mu(k,i)*gUb2(l,j)
 +#endif
              enddo
            enddo  
  cd         write (iout,*) 'EELEC: i',i,' j',j
@@@ -3366,55 -3252,10 +3380,55 @@@ cgrad            endi
  C Contribution to the local-electrostatic energy coming from the i-j pair
            eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
       &     +a33*muij(4)
 +c          write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4)
 +C Calculate patrial derivative for theta angle
 +#ifdef NEWCORR
 +         geel_loc_ij=a22*gmuij1(1)
 +     &     +a23*gmuij1(2)
 +     &     +a32*gmuij1(3)
 +     &     +a33*gmuij1(4)         
 +c         write(iout,*) "derivative over thatai"
 +c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
 +c     &   a33*gmuij1(4) 
 +         gloc(nphi+i,icg)=gloc(nphi+i,icg)+
 +     &      geel_loc_ij*wel_loc
 +c         write(iout,*) "derivative over thatai-1" 
 +c         write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
 +c     &   a33*gmuij2(4)
 +         geel_loc_ij=
 +     &     a22*gmuij2(1)
 +     &     +a23*gmuij2(2)
 +     &     +a32*gmuij2(3)
 +     &     +a33*gmuij2(4)
 +         gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
 +     &      geel_loc_ij*wel_loc
 +c  Derivative over j residue
 +         geel_loc_ji=a22*gmuji1(1)
 +     &     +a23*gmuji1(2)
 +     &     +a32*gmuji1(3)
 +     &     +a33*gmuji1(4)
 +c         write(iout,*) "derivative over thataj" 
 +c         write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
 +c     &   a33*gmuji1(4)
 +
 +        gloc(nphi+j,icg)=gloc(nphi+j,icg)+
 +     &      geel_loc_ji*wel_loc
 +         geel_loc_ji=
 +     &     +a22*gmuji2(1)
 +     &     +a23*gmuji2(2)
 +     &     +a32*gmuji2(3)
 +     &     +a33*gmuji2(4)
 +c         write(iout,*) "derivative over thataj-1"
 +c         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
 +c     &   a33*gmuji2(4)
 +         gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
 +     &      geel_loc_ji*wel_loc
 +#endif
  cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
  
            if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
       &            'eelloc',i,j,eel_loc_ij
 +c              write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
  
            eel_loc=eel_loc+eel_loc_ij
  C Partial derivatives in virtual-bond dihedral angles gamma
@@@ -3662,9 -3503,7 +3676,9 @@@ C Third- and fourth-order contribution
        dimension ggg(3)
        double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
       &  e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
 -     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2)
 +     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),gpizda1(2,2),
 +     &  gpizda2(2,2),auxgmat1(2,2),auxgmatt1(2,2),
 +     &  auxgmat2(2,2),auxgmatt2(2,2)
        double precision agg(3,4),aggi(3,4),aggi1(3,4),
       &    aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
        common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
@@@ -3688,24 -3527,9 +3702,24 @@@
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   
  cd        call checkint_turn3(i,a_temp,eello_turn3_num)
          call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
 +c auxalary matices for theta gradient
 +c auxalary matrix for i+1 and constant i+2
 +        call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1))
 +c auxalary matrix for i+2 and constant i+1
 +        call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1))
          call transpose2(auxmat(1,1),auxmat1(1,1))
 +        call transpose2(auxgmat1(1,1),auxgmatt1(1,1))
 +        call transpose2(auxgmat2(1,1),auxgmatt2(1,1))
          call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
 +        call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
 +        call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
          eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
 +C Derivatives in theta
 +        gloc(nphi+i,icg)=gloc(nphi+i,icg)
 +     &  +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
 +        gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
 +     &  +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
 +
          if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
       &          'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
  cd        write (2,*) 'i,',i,' j',j,'eello_turn3',
@@@ -3779,11 -3603,7 +3793,11 @@@ C Third- and fourth-order contribution
        dimension ggg(3)
        double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
       &  e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
 -     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2)
 +     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2),
 +     &  auxgEvec1(2),auxgEvec2(2),auxgEvec3(2),
 +     &  gte1t(2,2),gte2t(2,2),gte3t(2,2),
 +     &  gte1a(2,2),gtae3(2,2),gtae3e2(2,2), ae3gte2(2,2),
 +     &  gtEpizda1(2,2),gtEpizda2(2,2),gtEpizda3(2,2)
        double precision agg(3,4),aggi(3,4),aggi1(3,4),
       &    aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
        common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
@@@ -3803,7 -3623,6 +3817,7 @@@
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   
  cd        call checkint_turn4(i,a_temp,eello_turn4_num)
  c        write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
 +c        write(iout,*)"WCHODZE W PROGRAM"
          a_temp(1,1)=a22
          a_temp(1,2)=a23
          a_temp(2,1)=a32
@@@ -3815,95 -3634,32 +3829,95 @@@ c        write(iout,*) "iti1",iti1," it
          call transpose2(EUg(1,1,i+1),e1t(1,1))
          call transpose2(Eug(1,1,i+2),e2t(1,1))
          call transpose2(Eug(1,1,i+3),e3t(1,1))
 +C Ematrix derivative in theta
 +        call transpose2(gtEUg(1,1,i+1),gte1t(1,1))
 +        call transpose2(gtEug(1,1,i+2),gte2t(1,1))
 +        call transpose2(gtEug(1,1,i+3),gte3t(1,1))
          call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
 +c       eta1 in derivative theta
 +        call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1))
          call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
 -        s1=scalar2(b1(1,iti2),auxvec(1))
 +c       auxgvec is derivative of Ub2 so i+3 theta
 +        call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1)) 
 +c       auxalary matrix of E i+1
 +        call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1))
 +c        s1=0.0
 +c        gs1=0.0    
 +        s1=scalar2(b1(1,i+2),auxvec(1))
 +c derivative of theta i+2 with constant i+3
 +        gs23=scalar2(gtb1(1,i+2),auxvec(1))
 +c derivative of theta i+2 with constant i+2
 +        gs32=scalar2(b1(1,i+2),auxgvec(1))
 +c derivative of E matix in theta of i+1
 +        gsE13=scalar2(b1(1,i+2),auxgEvec1(1))
 +
          call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
 +c       ea31 in derivative theta
 +        call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1))
          call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
 -        s2=scalar2(b1(1,iti1),auxvec(1))
 +c auxilary matrix auxgvec of Ub2 with constant E matirx
 +        call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1))
 +c auxilary matrix auxgEvec1 of E matix with Ub2 constant
 +        call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1))
 +
 +c        s2=0.0
 +c        gs2=0.0
 +        s2=scalar2(b1(1,i+1),auxvec(1))
 +c derivative of theta i+1 with constant i+3
 +        gs13=scalar2(gtb1(1,i+1),auxvec(1))
 +c derivative of theta i+2 with constant i+1
 +        gs21=scalar2(b1(1,i+1),auxgvec(1))
 +c derivative of theta i+3 with constant i+1
 +        gsE31=scalar2(b1(1,i+1),auxgEvec3(1))
 +c        write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),gtb1(1,i+2),
 +c     &  gtb1(1,i+1)
          call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
 +c two derivatives over diffetent matrices
 +c gtae3e2 is derivative over i+3
 +        call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1))
 +c ae3gte2 is derivative over i+2
 +        call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1))
          call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
 +c three possible derivative over theta E matices
 +c i+1
 +        call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1))
 +c i+2
 +        call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1))
 +c i+3
 +        call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1))
          s3=0.5d0*(pizda(1,1)+pizda(2,2))
 +
 +        gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2))
 +        gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2))
 +        gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2))
 +
          eello_turn4=eello_turn4-(s1+s2+s3)
 +#ifdef NEWCORR
 +        gloc(nphi+i,icg)=gloc(nphi+i,icg)
 +     &                  -(gs13+gsE13+gsEE1)*wturn4
 +        gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)
 +     &                    -(gs23+gs21+gsEE2)*wturn4
 +        gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)
 +     &                    -(gs32+gsE31+gsEE3)*wturn4
 +c         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
 +c     &   gs2
 +#endif
          if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
       &      'eturn4',i,j,-(s1+s2+s3)
 -cd        write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
 -cd     &    ' eello_turn4_num',8*eello_turn4_num
 +c        write (iout,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
 +c     &    ' eello_turn4_num',8*eello_turn4_num
  C Derivatives in gamma(i)
          call transpose2(EUgder(1,1,i+1),e1tder(1,1))
          call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
          call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
 -        s1=scalar2(b1(1,iti2),auxvec(1))
 +        s1=scalar2(b1(1,i+2),auxvec(1))
          call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1))
          s3=0.5d0*(pizda(1,1)+pizda(2,2))
          gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3)
  C Derivatives in gamma(i+1)
          call transpose2(EUgder(1,1,i+2),e2tder(1,1))
          call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) 
 -        s2=scalar2(b1(1,iti1),auxvec(1))
 +        s2=scalar2(b1(1,i+1),auxvec(1))
          call matmat2(ae3(1,1),e2tder(1,1),auxmat(1,1))
          call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
          s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@@ -3911,10 -3667,10 +3925,10 @@@
  C Derivatives in gamma(i+2)
          call transpose2(EUgder(1,1,i+3),e3tder(1,1))
          call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1))
 -        s1=scalar2(b1(1,iti2),auxvec(1))
 +        s1=scalar2(b1(1,i+2),auxvec(1))
          call matmat2(a_temp(1,1),e3tder(1,1),auxmat(1,1))
          call matvec2(auxmat(1,1),Ub2(1,i+2),auxvec(1)) 
 -        s2=scalar2(b1(1,iti1),auxvec(1))
 +        s2=scalar2(b1(1,i+1),auxvec(1))
          call matmat2(auxmat(1,1),e2t(1,1),auxmat3(1,1))
          call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1))
          s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@@ -3929,10 -3685,10 +3943,10 @@@ C Derivatives of this turn contribution
              a_temp(2,2)=agg(l,4)
              call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
              call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
 -            s1=scalar2(b1(1,iti2),auxvec(1))
 +            s1=scalar2(b1(1,i+2),auxvec(1))
              call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
              call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
 -            s2=scalar2(b1(1,iti1),auxvec(1))
 +            s2=scalar2(b1(1,i+1),auxvec(1))
              call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
              call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
              s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@@ -3948,10 -3704,10 +3962,10 @@@ C Remaining derivatives of this turn co
            a_temp(2,2)=aggi(l,4)
            call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
            call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
 -          s1=scalar2(b1(1,iti2),auxvec(1))
 +          s1=scalar2(b1(1,i+2),auxvec(1))
            call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
            call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
 -          s2=scalar2(b1(1,iti1),auxvec(1))
 +          s2=scalar2(b1(1,i+1),auxvec(1))
            call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
            call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
            s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@@ -3962,10 -3718,10 +3976,10 @@@
            a_temp(2,2)=aggi1(l,4)
            call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
            call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
 -          s1=scalar2(b1(1,iti2),auxvec(1))
 +          s1=scalar2(b1(1,i+2),auxvec(1))
            call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
            call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
 -          s2=scalar2(b1(1,iti1),auxvec(1))
 +          s2=scalar2(b1(1,i+1),auxvec(1))
            call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
            call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
            s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@@ -3976,10 -3732,10 +3990,10 @@@
            a_temp(2,2)=aggj(l,4)
            call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
            call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
 -          s1=scalar2(b1(1,iti2),auxvec(1))
 +          s1=scalar2(b1(1,i+2),auxvec(1))
            call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
            call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
 -          s2=scalar2(b1(1,iti1),auxvec(1))
 +          s2=scalar2(b1(1,i+1),auxvec(1))
            call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
            call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
            s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@@ -3990,10 -3746,10 +4004,10 @@@
            a_temp(2,2)=aggj1(l,4)
            call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
            call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
 -          s1=scalar2(b1(1,iti2),auxvec(1))
 +          s1=scalar2(b1(1,i+2),auxvec(1))
            call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
            call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
 -          s2=scalar2(b1(1,iti1),auxvec(1))
 +          s2=scalar2(b1(1,i+1),auxvec(1))
            call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
            call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
            s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@@ -4060,7 -3816,7 +4074,7 @@@
  cd    print '(a)','Enter ESCP'
  cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
 -        if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
 +        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
          iteli=itel(i)
          xi=0.5D0*(c(1,i)+c(1,i+1))
          yi=0.5D0*(c(2,i)+c(2,i+1))
@@@ -4069,8 -3825,8 +4083,8 @@@
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
 -          if (itype(j).eq.21) cycle
 -          itypj=itype(j)
 +          if (itype(j).eq.ntyp1) cycle
 +          itypj=iabs(itype(j))
  C Uncomment following three lines for SC-p interactions
  c         xj=c(1,nres+j)-xi
  c         yj=c(2,nres+j)-yi
@@@ -4156,7 -3912,7 +4170,7 @@@
  cd    print '(a)','Enter ESCP'
  cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
        do i=iatscp_s,iatscp_e
 -        if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
 +        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
          iteli=itel(i)
          xi=0.5D0*(c(1,i)+c(1,i+1))
          yi=0.5D0*(c(2,i)+c(2,i+1))
@@@ -4165,8 -3921,8 +4179,8 @@@
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
 -          itypj=itype(j)
 -          if (itypj.eq.21) cycle
 +          itypj=iabs(itype(j))
 +          if (itypj.eq.ntyp1) cycle
  C Uncomment following three lines for SC-p interactions
  c         xj=c(1,nres+j)-xi
  c         yj=c(2,nres+j)-yi
@@@ -4186,9 -3942,8 +4200,9 @@@ C Uncomment following three lines for C
            endif
            evdwij=e1+e2
            evdw2=evdw2+evdwij
 -          if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 -     &        'evdw2',i,j,evdwij
 +          if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
 +     &        'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
 +     &       bad(itypj,iteli)
  C
  C Calculate contributions to the gradient in the virtual-bond and SC vectors.
  C
@@@ -4280,50 -4035,56 +4294,56 @@@ C iii and jjj point to the residues fo
            iii=ii
            jjj=jj
          endif
- cd        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
+ c        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
+ c     &    dhpb(i),dhpb1(i),forcon(i)
  C 24/11/03 AL: SS bridges handled separately because of introducing a specific
  C    distance and angle dependent SS bond potential.
 +        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
 +     & iabs(itype(jjj)).eq.1) then
+ cmc        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+ C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
+         if (.not.dyn_ss .and. i.le.nss) then
+ C 15/02/13 CC dynamic SSbond - additional check
 -         if (ii.gt.nres 
 -     &       .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then 
            call ssbond_ene(iii,jjj,eij)
            ehpb=ehpb+2*eij
+          endif
  cd          write (iout,*) "eij",eij
          else
  C Calculate the distance between the two points and its difference from the
  C target distance.
-         dd=dist(ii,jj)
-         rdis=dd-dhpb(i)
+           dd=dist(ii,jj)
+             rdis=dd-dhpb(i)
  C Get the force constant corresponding to this distance.
-         waga=forcon(i)
+             waga=forcon(i)
  C Calculate the contribution to energy.
-         ehpb=ehpb+waga*rdis*rdis
+             ehpb=ehpb+waga*rdis*rdis
  C
  C Evaluate gradient.
  C
-         fac=waga*rdis/dd
+             fac=waga*rdis/dd
  cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
  cd   &   ' waga=',waga,' fac=',fac
-         do j=1,3
-           ggg(j)=fac*(c(j,jj)-c(j,ii))
-         enddo
+             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
+           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
-         endif
+           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
+           do k=1,3
+             ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+             ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+           enddo
          endif
        enddo
        ehpb=0.5D0*ehpb
@@@ -4348,7 -4109,7 +4368,7 @@@
        include 'COMMON.VAR'
        include 'COMMON.IOUNITS'
        double precision erij(3),dcosom1(3),dcosom2(3),gg(3)
 -      itypi=itype(i)
 +      itypi=iabs(itype(i))
        xi=c(1,nres+i)
        yi=c(2,nres+i)
        zi=c(3,nres+i)
@@@ -4357,7 -4118,7 +4377,7 @@@
        dzi=dc_norm(3,nres+i)
  c      dsci_inv=dsc_inv(itypi)
        dsci_inv=vbld_inv(nres+i)
 -      itypj=itype(j)
 +      itypj=iabs(itype(j))
  c      dscj_inv=dsc_inv(itypj)
        dscj_inv=vbld_inv(nres+j)
        xj=c(1,nres+j)-xi
@@@ -4439,7 -4200,7 +4459,7 @@@
        estr=0.0d0
        estr1=0.0d0
        do i=ibondp_start,ibondp_end
 -        if (itype(i-1).eq.21 .or. itype(i).eq.21) then
 +        if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
            estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
            do j=1,3
            gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
@@@ -4463,8 -4224,8 +4483,8 @@@
  c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
  c
        do i=ibond_start,ibond_end
 -        iti=itype(i)
 -        if (iti.ne.10 .and. iti.ne.21) then
 +        iti=iabs(itype(i))
 +        if (iti.ne.10 .and. iti.ne.ntyp1) then
            nbi=nbondterm(iti)
            if (nbi.eq.1) then
              diff=vbld(i+nres)-vbldsc0(1,iti)
@@@ -4537,24 -4298,11 +4557,24 @@@ c      time12=1.0d
        etheta=0.0D0
  c     write (*,'(a,i2)') 'EBEND ICG=',icg
        do i=ithet_start,ithet_end
 -        if (itype(i-1).eq.21) cycle
 +        if (itype(i-1).eq.ntyp1) cycle
  C Zero the energy function and its derivative at 0 or pi.
          call splinthet(theta(i),0.5d0*delta,ss,ssd)
          it=itype(i-1)
 -        if (i.gt.3 .and. itype(i-2).ne.21) then
 +        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))
 +         endif
 +
 +        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
  #ifdef OSF
  	  phii=phi(i)
            if (phii.ne.phii) phii=150.0
@@@ -4567,7 -4315,7 +4587,7 @@@
            y(1)=0.0D0
            y(2)=0.0D0
          endif
 -        if (i.lt.nres .and. itype(i).ne.21) then
 +        if (i.lt.nres .and. itype(i).ne.ntyp1) then
  #ifdef OSF
  	  phii1=phi(i+1)
            if (phii1.ne.phii1) phii1=150.0
@@@ -4587,27 -4335,15 +4607,27 @@@ C dependent on the adjacent virtual-bon
  C In following comments this theta will be referred to as t_c.
          thet_pred_mean=0.0d0
          do k=1,2
 -          athetk=athet(k,it)
 -          bthetk=bthet(k,it)
 -          thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
 +            athetk=athet(k,it,ichir1,ichir2)
 +            bthetk=bthet(k,it,ichir1,ichir2)
 +          if (it.eq.10) then
 +             athetk=athet(k,itype1,ichir11,ichir12)
 +             bthetk=bthet(k,itype2,ichir21,ichir22)
 +          endif
 +         thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
          enddo
          dthett=thet_pred_mean*ssd
          thet_pred_mean=thet_pred_mean*ss+a0thet(it)
  C Derivatives of the "mean" values in gamma1 and gamma2.
 -        dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss
 -        dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss
 +        dthetg1=(-athet(1,it,ichir1,ichir2)*y(2)
 +     &+athet(2,it,ichir1,ichir2)*y(1))*ss
 +         dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2)
 +     &          +bthet(2,it,ichir1,ichir2)*z(1))*ss
 +         if (it.eq.10) then
 +      dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2)
 +     &+athet(2,itype1,ichir11,ichir12)*y(1))*ss
 +        dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2)
 +     &         +bthet(2,itype2,ichir21,ichir22)*z(1))*ss
 +         endif
          if (theta(i).gt.pi-delta) then
            call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
       &         E_tc0)
@@@ -4634,7 -4370,7 +4654,7 @@@
       &      'ebend',i,ethetai
          if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
          if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
 -        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
 +        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
        enddo
  C Ufff.... We've done all this!!! 
        return
@@@ -4775,27 -4511,24 +4795,27 @@@
        logical lprn /.false./, lprn1 /.false./
        etheta=0.0D0
        do i=ithet_start,ithet_end
 -        if (itype(i-1).eq.21) cycle
 +        if (itype(i-1).eq.ntyp1) cycle
 +        if (iabs(itype(i+1)).eq.20) iblock=2
 +        if (iabs(itype(i+1)).ne.20) iblock=1
          dethetai=0.0d0
          dephii=0.0d0
          dephii1=0.0d0
          theti2=0.5d0*theta(i)
 -        ityp2=ithetyp(itype(i-1))
 +        ityp2=ithetyp((itype(i-1)))
          do k=1,nntheterm
            coskt(k)=dcos(k*theti2)
            sinkt(k)=dsin(k*theti2)
          enddo
 -        if (i.gt.3 .and. itype(i-2).ne.21) then
 +        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
  #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)))
 +C propagation of chirality for glycine type
            do k=1,nsingle
              cosph1(k)=dcos(k*phii)
              sinph1(k)=dsin(k*phii)
@@@ -4808,7 -4541,7 +4828,7 @@@
              sinph1(k)=0.0d0
            enddo 
          endif
 -        if (i.lt.nres .and. itype(i).ne.21) then
 +        if (i.lt.nres .and. itype(i).ne.ntyp1) then
  #ifdef OSF
            phii1=phi(i+1)
            if (phii1.ne.phii1) phii1=150.0
@@@ -4816,7 -4549,7 +4836,7 @@@
  #else
            phii1=phi(i+1)
  #endif
 -          ityp3=ithetyp(itype(i))
 +          ityp3=ithetyp((itype(i)))
            do k=1,nsingle
              cosph2(k)=dcos(k*phii1)
              sinph2(k)=dsin(k*phii1)
@@@ -4829,7 -4562,7 +4849,7 @@@
              sinph2(k)=0.0d0
            enddo
          endif  
 -        ethetai=aa0thet(ityp1,ityp2,ityp3)
 +        ethetai=aa0thet(ityp1,ityp2,ityp3,iblock)
          do k=1,ndouble
            do l=1,k-1
              ccl=cosph1(l)*cosph2(k-l)
@@@ -4851,12 -4584,11 +4871,12 @@@
          enddo
          endif
          do k=1,ntheterm
 -          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
 -          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
 +          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k)
 +          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock)
       &      *coskt(k)
            if (lprn)
 -     &    write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
 +     &    write (iout,*) "k",k,"
 +     &     aathet",aathet(k,ityp1,ityp2,ityp3,iblock),
       &     " ethetai",ethetai
          enddo
          if (lprn) then
@@@ -4875,24 -4607,24 +4895,24 @@@
          endif
          do m=1,ntheterm2
            do k=1,nsingle
 -            aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
 -     &         +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
 -     &         +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
 -     &         +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
 +            aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
 +     &         +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)
 +     &         +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)
 +     &         +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)
              ethetai=ethetai+sinkt(m)*aux
              dethetai=dethetai+0.5d0*m*aux*coskt(m)
              dephii=dephii+k*sinkt(m)*(
 -     &          ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
 -     &          bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
 +     &          ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)-
 +     &          bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k))
              dephii1=dephii1+k*sinkt(m)*(
 -     &          eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
 -     &          ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
 +     &          eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)-
 +     &          ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k))
              if (lprn)
       &      write (iout,*) "m",m," k",k," bbthet",
 -     &         bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
 -     &         ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
 -     &         ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
 -     &         eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
 +     &         bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet",
 +     &         ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
 +     &         ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
 +     &         eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
            enddo
          enddo
          if (lprn)
@@@ -4900,29 -4632,28 +4920,29 @@@
          do m=1,ntheterm3
            do k=2,ndouble
              do l=1,k-1
 -              aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
 +              aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
 +     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+
 +     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
 +     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)
                ethetai=ethetai+sinkt(m)*aux
                dethetai=dethetai+0.5d0*m*coskt(m)*aux
                dephii=dephii+l*sinkt(m)*(
 -     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
 +     &           -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)-
 +     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
 +     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
 +     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
                dephii1=dephii1+(k-l)*sinkt(m)*(
 -     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
 +     &           -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
 +     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
 +     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)-
 +     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
                if (lprn) then
                write (iout,*) "m",m," k",k," l",l," ffthet",
 -     &            ffthet(l,k,m,ityp1,ityp2,ityp3),
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3),
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
 +     &            ffthet(l,k,m,ityp1,ityp2,ityp3,iblock),
 +     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet",
 +     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock),
 +     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock),
 +     &            " ethetai",ethetai
                write (iout,*) cosph1ph2(l,k)*sinkt(m),
       &            cosph1ph2(k,l)*sinkt(m),
       &            sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
@@@ -4931,16 -4662,13 +4951,16 @@@
            enddo
          enddo
  10      continue
 -        if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
 +c        lprn1=.true.
 +        if (lprn1) 
 +     &    write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
       &   i,theta(i)*rad2deg,phii*rad2deg,
       &   phii1*rad2deg,ethetai
 +c        lprn1=.false.
          etheta=etheta+ethetai
          if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
          if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
 -        gloc(nphi+i-2,icg)=wang*dethetai
 +        gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg)
        enddo
        return
        end
@@@ -4971,9 -4699,9 +4991,9 @@@ C ALPHA and OMEGA
  c     write (iout,'(a)') 'ESC'
        do i=loc_start,loc_end
          it=itype(i)
 -        if (it.eq.21) cycle
 +        if (it.eq.ntyp1) cycle
          if (it.eq.10) goto 1
 -        nlobit=nlob(it)
 +        nlobit=nlob(iabs(it))
  c       print *,'i=',i,' it=',it,' nlobit=',nlobit
  c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
          theti=theta(i+1)-pipol
@@@ -5130,11 -4858,11 +5150,11 @@@ C Compute the contribution to SC energ
  
            do j=1,nlobit
  #ifdef OSF
 -            adexp=bsc(j,it)-0.5D0*contr(j,iii)+emin
 +            adexp=bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin
              if(adexp.ne.adexp) adexp=1.0
              expfac=dexp(adexp)
  #else
 -            expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin)
 +            expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin)
  #endif
  cd          print *,'j=',j,' expfac=',expfac
              escloc_i=escloc_i+expfac
@@@ -5216,7 -4944,7 +5236,7 @@@ C Compute the contribution to SC energ
  
        dersc12=0.0d0
        do j=1,nlobit
 -        expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin)
 +        expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin)
          escloc_i=escloc_i+expfac
          do k=1,2
            dersc(k)=dersc(k)+Ax(k,j)*expfac
@@@ -5270,7 -4998,7 +5290,7 @@@
        delta=0.02d0*pi
        escloc=0.0D0
        do i=loc_start,loc_end
 -        if (itype(i).eq.21) cycle
 +        if (itype(i).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)))
@@@ -5279,7 -5007,7 +5299,7 @@@
          cosfac=dsqrt(cosfac2)
          sinfac2=0.5d0/(1.0d0-costtab(i+1))
          sinfac=dsqrt(sinfac2)
 -        it=itype(i)
 +        it=iabs(itype(i))
          if (it.eq.10) goto 1
  c
  C  Compute the axes of tghe local cartesian coordinates system; store in
@@@ -5297,7 -5025,7 +5317,7 @@@ C     &   dc_norm(3,i+nres
            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)
 +          z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
          enddo     
  c       write (2,*) "i",i
  c       write (2,*) "x_prime",(x_prime(j),j=1,3)
@@@ -5329,7 -5057,7 +5349,7 @@@
  C Compute the energy of the ith side cbain
  C
  c        write (2,*) "xx",xx," yy",yy," zz",zz
 -        it=itype(i)
 +        it=iabs(itype(i))
          do j = 1,65
            x(j) = sc_parmin(j,it) 
          enddo
@@@ -5337,7 -5065,7 +5357,7 @@@
  Cc diagnostics - remove later
          xx1 = dcos(alph(2))
          yy1 = dsin(alph(2))*dcos(omeg(2))
 -        zz1 = -dsin(alph(2))*dsin(omeg(2))
 +        zz1 = -dsign(1.0,dfloat(itype(i)))*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
@@@ -5379,9 -5107,7 +5399,9 @@@ c     &   sumene4
  c     &   dscp1,dscp2,sumene
  c        sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
          escloc = escloc + sumene
 -c        write (2,*) "i",i," escloc",sumene,escloc
 +c        write (2,*) "i",i," escloc",sumene,escloc,it,itype(i)
 +c     & ,zz,xx,yy
 +c#define DEBUG
  #ifdef DEBUG
  C
  C This section to check the numerical derivatives of the energy of ith side
@@@ -5425,7 -5151,6 +5445,7 @@@ C End of diagnostics section
  C        
  C Compute the gradient of esc
  C
 +c        zz=zz*dsign(1.0,dfloat(itype(i)))
          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
@@@ -5450,7 -5175,7 +5470,7 @@@
       &        +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6)
       &        +(pom1+pom2)*pom_dx
  #ifdef DEBUG
 -        write(2,*), "de_dxx = ", de_dxx,de_dxx_num
 +        write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i)
  #endif
  C
          sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz
@@@ -5465,7 -5190,7 +5485,7 @@@
       &        +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6)
       &        +(pom1-pom2)*pom_dy
  #ifdef DEBUG
 -        write(2,*), "de_dyy = ", de_dyy,de_dyy_num
 +        write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i)
  #endif
  C
          de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy
@@@ -5477,16 -5202,15 +5497,16 @@@
       &  +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
 +        write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i)
  #endif
  C
          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
 +        write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i)
  #endif
 +c#undef DEBUG
  c 
  C
         cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
@@@ -5511,16 -5235,13 +5531,16 @@@ c     &   (dC_norm(j,i-1),j=1,3)," vbld
           dZZ_Ci1(k)=0.0d0
           dZZ_Ci(k)=0.0d0
           do j=1,3
 -           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
 -           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
 +           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
 +     &     *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
 +           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)
 +     &     *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
           enddo
            
           dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
           dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres))
 -         dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres))
 +         dZZ_XYZ(k)=vbld_inv(i+nres)*
 +     &   (z_prime(k)-zz*dC_norm(k,i+nres))
  c
           dt_dCi(k) = -dt_dCi(k)/sinttab(i+1)
           dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1)
@@@ -5705,8 -5426,8 +5725,8 @@@ c      lprn=.true
        etors=0.0D0
        do i=iphi_start,iphi_end
        etors_ii=0.0D0
 -        if (itype(i-2).eq.21 .or. itype(i-1).eq.21 
 -     &      .or. itype(i).eq.21) cycle
 +        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))
          phii=phi(i)
@@@ -5802,22 -5523,17 +5822,22 @@@ C Set lprn=.true. for debuggin
  c     lprn=.true.
        etors=0.0D0
        do i=iphi_start,iphi_end
 -        if (itype(i-2).eq.21 .or. itype(i-1).eq.21 
 -     &       .or. itype(i).eq.21) cycle
 -      etors_ii=0.0D0
 +        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 
 +     &       .or. itype(i).eq.ntyp1) cycle
 +        etors_ii=0.0D0
 +         if (iabs(itype(i)).eq.20) then
 +         iblock=2
 +         else
 +         iblock=1
 +         endif
          itori=itortyp(itype(i-2))
          itori1=itortyp(itype(i-1))
          phii=phi(i)
          gloci=0.0D0
  C Regular cosine and sine terms
 -        do j=1,nterm(itori,itori1)
 -          v1ij=v1(j,itori,itori1)
 -          v2ij=v2(j,itori,itori1)
 +        do j=1,nterm(itori,itori1,iblock)
 +          v1ij=v1(j,itori,itori1,iblock)
 +          v2ij=v2(j,itori,itori1,iblock)
            cosphi=dcos(j*phii)
            sinphi=dsin(j*phii)
            etors=etors+v1ij*cosphi+v2ij*sinphi
@@@ -5832,7 -5548,7 +5852,7 @@@ C          [v2 cos(phi/2)+v3 sin(phi/2)
  C
          cosphi=dcos(0.5d0*phii)
          sinphi=dsin(0.5d0*phii)
 -        do j=1,nlor(itori,itori1)
 +        do j=1,nlor(itori,itori1,iblock)
            vl1ij=vlor1(j,itori,itori1)
            vl2ij=vlor2(j,itori,itori1)
            vl3ij=vlor3(j,itori,itori1)
@@@ -5845,14 -5561,13 +5865,14 @@@
            gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
          enddo
  C Subtract the constant term
 -        etors=etors-v0(itori,itori1)
 +        etors=etors-v0(itori,itori1,iblock)
            if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
 -     &         'etor',i,etors_ii-v0(itori,itori1)
 +     &         '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,
 -     &  (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
 +     &  (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
  c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
        enddo
@@@ -5902,10 -5617,9 +5922,10 @@@ C Set lprn=.true. for debuggin
        lprn=.false.
  c     lprn=.true.
        etors_d=0.0D0
 +c      write(iout,*) "a tu??"
        do i=iphid_start,iphid_end
 -        if (itype(i-2).eq.21 .or. itype(i-1).eq.21
 -     &      .or. itype(i).eq.21 .or. itype(i+1).eq.21) cycle
 +        if (itype(i-2).eq.ntyp1 .or. itype(i-1).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))
@@@ -5913,15 -5627,12 +5933,15 @@@
          phii1=phi(i+1)
          gloci1=0.0D0
          gloci2=0.0D0
 +        iblock=1
 +        if (iabs(itype(i+1)).eq.20) iblock=2
 +
  C Regular cosine and sine terms
 -        do j=1,ntermd_1(itori,itori1,itori2)
 -          v1cij=v1c(1,j,itori,itori1,itori2)
 -          v1sij=v1s(1,j,itori,itori1,itori2)
 -          v2cij=v1c(2,j,itori,itori1,itori2)
 -          v2sij=v1s(2,j,itori,itori1,itori2)
 +        do j=1,ntermd_1(itori,itori1,itori2,iblock)
 +          v1cij=v1c(1,j,itori,itori1,itori2,iblock)
 +          v1sij=v1s(1,j,itori,itori1,itori2,iblock)
 +          v2cij=v1c(2,j,itori,itori1,itori2,iblock)
 +          v2sij=v1s(2,j,itori,itori1,itori2,iblock)
            cosphi1=dcos(j*phii)
            sinphi1=dsin(j*phii)
            cosphi2=dcos(j*phii1)
@@@ -5931,12 -5642,12 +5951,12 @@@
            gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
            gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
          enddo
 -        do k=2,ntermd_2(itori,itori1,itori2)
 +        do k=2,ntermd_2(itori,itori1,itori2,iblock)
            do l=1,k-1
 -            v1cdij = v2c(k,l,itori,itori1,itori2)
 -            v2cdij = v2c(l,k,itori,itori1,itori2)
 -            v1sdij = v2s(k,l,itori,itori1,itori2)
 -            v2sdij = v2s(l,k,itori,itori1,itori2)
 +            v1cdij = v2c(k,l,itori,itori1,itori2,iblock)
 +            v2cdij = v2c(l,k,itori,itori1,itori2,iblock)
 +            v1sdij = v2s(k,l,itori,itori1,itori2,iblock)
 +            v2sdij = v2s(l,k,itori,itori1,itori2,iblock)
              cosphi1p2=dcos(l*phii+(k-l)*phii1)
              cosphi1m2=dcos(l*phii-(k-l)*phii1)
              sinphi1p2=dsin(l*phii+(k-l)*phii1)
@@@ -5981,53 -5692,29 +6001,53 @@@ c        amino-acid residues
  C Set lprn=.true. for debugging
        lprn=.false.
  c      lprn=.true.
 -c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
 +c      write (iout,*) "EBACK_SC_COR",itau_start,itau_end
        esccor=0.0D0
 -      do i=iphi_start,iphi_end
 -        if (itype(i-2).eq.21 .or. itype(i-1).eq.21) cycle
 +      do i=itau_start,itau_end
 +        if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
          esccor_ii=0.0D0
 -        itori=itype(i-2)
 -        itori1=itype(i-1)
 +        isccori=isccortyp(itype(i-2))
 +        isccori1=isccortyp(itype(i-1))
 +c      write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
          phii=phi(i)
 +        do intertyp=1,3 !intertyp
 +cc Added 09 May 2012 (Adasko)
 +cc  Intertyp means interaction type of backbone mainchain correlation: 
 +c   1 = SC...Ca...Ca...Ca
 +c   2 = Ca...Ca...Ca...SC
 +c   3 = SC...Ca...Ca...SCi
          gloci=0.0D0
 -        do j=1,nterm_sccor
 -          v1ij=v1sccor(j,itori,itori1)
 -          v2ij=v2sccor(j,itori,itori1)
 -          cosphi=dcos(j*phii)
 -          sinphi=dsin(j*phii)
 +        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))
 +     & cycle
 +       do j=1,nterm_sccor(isccori,isccori1)
 +          v1ij=v1sccor(j,intertyp,isccori,isccori1)
 +          v2ij=v2sccor(j,intertyp,isccori,isccori1)
 +          cosphi=dcos(j*tauangle(intertyp,i))
 +          sinphi=dsin(j*tauangle(intertyp,i))
            esccor=esccor+v1ij*cosphi+v2ij*sinphi
            gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
          enddo
 +c      write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
 +        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,itori,itori1,
 -     &  (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
 +     &  restyp(itype(i-2)),i-2,restyp(itype(i-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
 +       enddo !intertyp
        enddo
 +
        return
        end
  c----------------------------------------------------------------------------
@@@ -7037,10 -6724,10 +7057,10 @@@ C--------------------------------------
        do iii=1,2
          dipi(iii,1)=Ub2(iii,i)
          dipderi(iii)=Ub2der(iii,i)
 -        dipi(iii,2)=b1(iii,iti1)
 +        dipi(iii,2)=b1(iii,i+1)
          dipj(iii,1)=Ub2(iii,j)
          dipderj(iii)=Ub2der(iii,j)
 -        dipj(iii,2)=b1(iii,itj1)
 +        dipj(iii,2)=b1(iii,j+1)
        enddo
        kkk=0
        do iii=1,2
@@@ -7220,26 -6907,26 +7240,26 @@@ C They are needed only when the fifth- 
  C indluded.
          IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) THEN
          call transpose2(AEA(1,1,1),auxmat(1,1))
 -        call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
 +        call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1))
          call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
          call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
          call transpose2(AEAderg(1,1,1),auxmat(1,1))
 -        call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
 +        call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1))
          call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
 -        call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
 -        call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
 +        call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1))
 +        call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1))
          call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
          call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
          call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
          call transpose2(AEA(1,1,2),auxmat(1,1))
 -        call matvec2(auxmat(1,1),b1(1,itj),AEAb1(1,1,2))
 +        call matvec2(auxmat(1,1),b1(1,j),AEAb1(1,1,2))
          call matvec2(auxmat(1,1),Ub2(1,j),AEAb2(1,1,2))
          call matvec2(auxmat(1,1),Ub2der(1,j),AEAb2derg(1,2,1,2))
          call transpose2(AEAderg(1,1,2),auxmat(1,1))
 -        call matvec2(auxmat(1,1),b1(1,itj),AEAb1derg(1,1,2))
 +        call matvec2(auxmat(1,1),b1(1,j),AEAb1derg(1,1,2))
          call matvec2(auxmat(1,1),Ub2(1,j),AEAb2derg(1,1,1,2))
 -        call matvec2(AEA(1,1,2),b1(1,itl1),AEAb1(1,2,2))
 -        call matvec2(AEAderg(1,1,2),b1(1,itl1),AEAb1derg(1,2,2))
 +        call matvec2(AEA(1,1,2),b1(1,l+1),AEAb1(1,2,2))
 +        call matvec2(AEAderg(1,1,2),b1(1,l+1),AEAb1derg(1,2,2))
          call matvec2(AEA(1,1,2),Ub2(1,l+1),AEAb2(1,2,2))
          call matvec2(AEAderg(1,1,2),Ub2(1,l+1),AEAb2derg(1,1,2,2))
          call matvec2(AEA(1,1,2),Ub2der(1,l+1),AEAb2derg(1,2,2,2))
@@@ -7248,20 -6935,20 +7268,20 @@@ C Calculate the Cartesian derivatives o
            do kkk=1,5
              do lll=1,3
                call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
 -              call matvec2(auxmat(1,1),b1(1,iti),
 +              call matvec2(auxmat(1,1),b1(1,i),
       &          AEAb1derx(1,lll,kkk,iii,1,1))
                call matvec2(auxmat(1,1),Ub2(1,i),
       &          AEAb2derx(1,lll,kkk,iii,1,1))
 -              call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
 +              call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
       &          AEAb1derx(1,lll,kkk,iii,2,1))
                call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
       &          AEAb2derx(1,lll,kkk,iii,2,1))
                call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
 -              call matvec2(auxmat(1,1),b1(1,itj),
 +              call matvec2(auxmat(1,1),b1(1,j),
       &          AEAb1derx(1,lll,kkk,iii,1,2))
                call matvec2(auxmat(1,1),Ub2(1,j),
       &          AEAb2derx(1,lll,kkk,iii,1,2))
 -              call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
 +              call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,l+1),
       &          AEAb1derx(1,lll,kkk,iii,2,2))
                call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,l+1),
       &          AEAb2derx(1,lll,kkk,iii,2,2))
@@@ -7358,26 -7045,26 +7378,26 @@@ C indluded
          IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 .or.
       &    (wturn6.gt.0.0d0 .and. j.eq.i+4 .and. l.eq.i+3)) THEN
          call transpose2(AEA(1,1,1),auxmat(1,1))
 -        call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
 +        call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1))
          call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
          call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
          call transpose2(AEAderg(1,1,1),auxmat(1,1))
 -        call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
 +        call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1))
          call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
 -        call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
 -        call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
 +        call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1))
 +        call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1))
          call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
          call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
          call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
          call transpose2(AEA(1,1,2),auxmat(1,1))
 -        call matvec2(auxmat(1,1),b1(1,itj1),AEAb1(1,1,2))
 +        call matvec2(auxmat(1,1),b1(1,j+1),AEAb1(1,1,2))
          call matvec2(auxmat(1,1),Ub2(1,l),AEAb2(1,1,2))
          call matvec2(auxmat(1,1),Ub2der(1,l),AEAb2derg(1,2,1,2))
          call transpose2(AEAderg(1,1,2),auxmat(1,1))
 -        call matvec2(auxmat(1,1),b1(1,itl),AEAb1(1,1,2))
 +        call matvec2(auxmat(1,1),b1(1,l),AEAb1(1,1,2))
          call matvec2(auxmat(1,1),Ub2(1,l),AEAb2derg(1,1,1,2))
 -        call matvec2(AEA(1,1,2),b1(1,itj1),AEAb1(1,2,2))
 -        call matvec2(AEAderg(1,1,2),b1(1,itj1),AEAb1derg(1,2,2))
 +        call matvec2(AEA(1,1,2),b1(1,j+1),AEAb1(1,2,2))
 +        call matvec2(AEAderg(1,1,2),b1(1,j+1),AEAb1derg(1,2,2))
          call matvec2(AEA(1,1,2),Ub2(1,j),AEAb2(1,2,2))
          call matvec2(AEAderg(1,1,2),Ub2(1,j),AEAb2derg(1,1,2,2))
          call matvec2(AEA(1,1,2),Ub2der(1,j),AEAb2derg(1,2,2,2))
@@@ -7386,20 -7073,20 +7406,20 @@@ C Calculate the Cartesian derivatives o
            do kkk=1,5
              do lll=1,3
                call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
 -              call matvec2(auxmat(1,1),b1(1,iti),
 +              call matvec2(auxmat(1,1),b1(1,i),
       &          AEAb1derx(1,lll,kkk,iii,1,1))
                call matvec2(auxmat(1,1),Ub2(1,i),
       &          AEAb2derx(1,lll,kkk,iii,1,1))
 -              call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
 +              call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
       &          AEAb1derx(1,lll,kkk,iii,2,1))
                call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
       &          AEAb2derx(1,lll,kkk,iii,2,1))
                call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
 -              call matvec2(auxmat(1,1),b1(1,itl),
 +              call matvec2(auxmat(1,1),b1(1,l),
       &          AEAb1derx(1,lll,kkk,iii,1,2))
                call matvec2(auxmat(1,1),Ub2(1,l),
       &          AEAb2derx(1,lll,kkk,iii,1,2))
 -              call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itj1),
 +              call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,j+1),
       &          AEAb1derx(1,lll,kkk,iii,2,2))
                call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,j),
       &          AEAb2derx(1,lll,kkk,iii,2,2))
@@@ -7696,7 -7383,7 +7716,7 @@@ C Contribution from graph I
        call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1))
        vv(1)=pizda(1,1)+pizda(2,2)
        vv(2)=pizda(2,1)-pizda(1,2)
 -      eello5_2=scalar2(AEAb1(1,2,1),b1(1,itk))
 +      eello5_2=scalar2(AEAb1(1,2,1),b1(1,k))
       & -0.5d0*scalar2(vv(1),Ctobr(1,k))
  C Explicit gradient in virtual-dihedral angles.
        g_corr5_loc(k-1)=g_corr5_loc(k-1)
@@@ -7706,11 -7393,11 +7726,11 @@@
        vv(2)=pizda(2,1)-pizda(1,2)
        if (l.eq.j+1) then
          g_corr5_loc(l-1)=g_corr5_loc(l-1)
 -     &   +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
 +     &   +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k))
       &   -0.5d0*scalar2(vv(1),Ctobr(1,k)))
        else
          g_corr5_loc(j-1)=g_corr5_loc(j-1)
 -     &   +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
 +     &   +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k))
       &   -0.5d0*scalar2(vv(1),Ctobr(1,k)))
        endif
  C Cartesian gradient
@@@ -7722,7 -7409,7 +7742,7 @@@
              vv(1)=pizda(1,1)+pizda(2,2)
              vv(2)=pizda(2,1)-pizda(1,2)
              derx(lll,kkk,iii)=derx(lll,kkk,iii)
 -     &       +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,itk))
 +     &       +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,k))
       &       -0.5d0*scalar2(vv(1),Ctobr(1,k))
            enddo
          enddo
@@@ -7777,7 -7464,7 +7797,7 @@@ cd1110    continu
          call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
          vv(1)=pizda(1,1)+pizda(2,2)
          vv(2)=pizda(2,1)-pizda(1,2)
 -        eello5_4=scalar2(AEAb1(1,2,2),b1(1,itl))
 +        eello5_4=scalar2(AEAb1(1,2,2),b1(1,l))
       &   -0.5d0*scalar2(vv(1),Ctobr(1,l))
  C Explicit gradient in virtual-dihedral angles.
          g_corr5_loc(l-1)=g_corr5_loc(l-1)
@@@ -7786,7 -7473,7 +7806,7 @@@
          vv(1)=pizda(1,1)+pizda(2,2)
          vv(2)=pizda(2,1)-pizda(1,2)
          g_corr5_loc(k-1)=g_corr5_loc(k-1)
 -     &   +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itl))
 +     &   +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,l))
       &   -0.5d0*scalar2(vv(1),Ctobr(1,l)))
  C Cartesian gradient
          do iii=1,2
@@@ -7797,7 -7484,7 +7817,7 @@@
                vv(1)=pizda(1,1)+pizda(2,2)
                vv(2)=pizda(2,1)-pizda(1,2)
                derx(lll,kkk,iii)=derx(lll,kkk,iii)
 -     &         +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itl))
 +     &         +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,l))
       &         -0.5d0*scalar2(vv(1),Ctobr(1,l))
              enddo
            enddo
@@@ -7850,7 -7537,7 +7870,7 @@@ C Contribution from graph I
          call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
          vv(1)=pizda(1,1)+pizda(2,2)
          vv(2)=pizda(2,1)-pizda(1,2)
 -        eello5_4=scalar2(AEAb1(1,2,2),b1(1,itj))
 +        eello5_4=scalar2(AEAb1(1,2,2),b1(1,j))
       &   -0.5d0*scalar2(vv(1),Ctobr(1,j))
  C Explicit gradient in virtual-dihedral angles.
          g_corr5_loc(j-1)=g_corr5_loc(j-1)
@@@ -7859,7 -7546,7 +7879,7 @@@
          vv(1)=pizda(1,1)+pizda(2,2)
          vv(2)=pizda(2,1)-pizda(1,2)
          g_corr5_loc(k-1)=g_corr5_loc(k-1)
 -     &   +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itj))
 +     &   +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,j))
       &   -0.5d0*scalar2(vv(1),Ctobr(1,j)))
  C Cartesian gradient
          do iii=1,2
@@@ -7870,7 -7557,7 +7890,7 @@@
                vv(1)=pizda(1,1)+pizda(2,2)
                vv(2)=pizda(2,1)-pizda(1,2)
                derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)
 -     &         +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itj))
 +     &         +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,j))
       &         -0.5d0*scalar2(vv(1),Ctobr(1,j))
              enddo
            enddo
@@@ -8152,8 -7839,8 +8172,8 @@@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        vv1(1)=pizda1(1,1)-pizda1(2,2)
        vv1(2)=pizda1(1,2)+pizda1(2,1)
        s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
 -      vv(1)=AEAb1(1,2,imat)*b1(1,itk)-AEAb1(2,2,imat)*b1(2,itk)
 -      vv(2)=AEAb1(1,2,imat)*b1(2,itk)+AEAb1(2,2,imat)*b1(1,itk)
 +      vv(1)=AEAb1(1,2,imat)*b1(1,k)-AEAb1(2,2,imat)*b1(2,k)
 +      vv(2)=AEAb1(1,2,imat)*b1(2,k)+AEAb1(2,2,imat)*b1(1,k)
        s5=scalar2(vv(1),Dtobr2(1,i))
  cd      write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5
        eello6_graph1=-0.5d0*(s1+s2+s3+s4+s5)
@@@ -8166,8 -7853,8 +8186,8 @@@
        call matmat2(AEAderg(1,1,imat),auxmat(1,1),pizda1(1,1))
        vv1(1)=pizda1(1,1)-pizda1(2,2)
        vv1(2)=pizda1(1,2)+pizda1(2,1)
 -      vv(1)=AEAb1derg(1,2,imat)*b1(1,itk)-AEAb1derg(2,2,imat)*b1(2,itk)
 -      vv(2)=AEAb1derg(1,2,imat)*b1(2,itk)+AEAb1derg(2,2,imat)*b1(1,itk)
 +      vv(1)=AEAb1derg(1,2,imat)*b1(1,k)-AEAb1derg(2,2,imat)*b1(2,k)
 +      vv(2)=AEAb1derg(1,2,imat)*b1(2,k)+AEAb1derg(2,2,imat)*b1(1,k)
        if (l.eq.j+1) then
          g_corr6_loc(l-1)=g_corr6_loc(l-1)
       & +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i))
@@@ -8206,10 -7893,10 +8226,10 @@@
              vv1(1)=pizda1(1,1)-pizda1(2,2)
              vv1(2)=pizda1(1,2)+pizda1(2,1)
              s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
 -            vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,itk)
 -     &       -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,itk)
 -            vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,itk)
 -     &       +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,itk)
 +            vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,k)
 +     &       -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,k)
 +            vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,k)
 +     &       +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,k)
              s5=scalar2(vv(1),Dtobr2(1,i))
              derx(lll,kkk,ind)=derx(lll,kkk,ind)-0.5d0*(s1+s2+s3+s4+s5)
            enddo
@@@ -8241,12 -7928,12 +8261,12 @@@
  C          o             o                                                     C
  C     \   /l\           /j\   /                                                C
  C      \ /   \         /   \ /                                                 C
- C       o| o |         | o |o                                                  C
+ C       o| o |         | o |o                                                  C                
  C     \ j|/k\|      \  |/k\|l                                                  C
  C      \ /   \       \ /   \                                                   C
  C       o             o                                                        C
- C       i             i                                                        C
- C                                                                              C
+ C       i             i                                                        C 
+ C                                                                              C           
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  cd      write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l
  C AL 7/4/01 s1 would occur in the sixth-order moment, 
@@@ -8417,10 -8104,10 +8437,10 @@@ c--------------------------------------
        double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
        logical swap
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- C                                                                              C
+ C                                                                              C 
  C      Parallel       Antiparallel                                             C
  C                                                                              C
- C          o             o                                                     C
+ C          o             o                                                     C 
  C         /l\   /   \   /j\                                                    C 
  C        /   \ /     \ /   \                                                   C
  C       /| o |o       o| o |\                                                  C
@@@ -8449,10 -8136,10 +8469,10 @@@ C           energy moment and not to th
  #ifdef MOMENT
        s1=dip(4,jj,i)*dip(4,kk,k)
  #endif
 -      call matvec2(AECA(1,1,1),b1(1,itk1),auxvec(1))
 -      s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
 -      call matvec2(AECA(1,1,2),b1(1,itl1),auxvec(1))
 -      s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
 +      call matvec2(AECA(1,1,1),b1(1,k+1),auxvec(1))
 +      s2=0.5d0*scalar2(b1(1,k),auxvec(1))
 +      call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1))
 +      s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
        call transpose2(EE(1,1,itk),auxmat(1,1))
        call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1))
        vv(1)=pizda(1,1)+pizda(2,2)
@@@ -8467,13 -8154,13 +8487,13 @@@ cd     & "sum",-(s2+s3+s4
  #endif
  c      eello6_graph3=-s4
  C Derivatives in gamma(k-1)
 -      call matvec2(AECAderg(1,1,2),b1(1,itl1),auxvec(1))
 -      s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
 +      call matvec2(AECAderg(1,1,2),b1(1,l+1),auxvec(1))
 +      s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
        s4=-0.25d0*scalar2(vv(1),Ctobrder(1,k))
        g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s3+s4)
  C Derivatives in gamma(l-1)
 -      call matvec2(AECAderg(1,1,1),b1(1,itk1),auxvec(1))
 -      s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
 +      call matvec2(AECAderg(1,1,1),b1(1,k+1),auxvec(1))
 +      s2=0.5d0*scalar2(b1(1,k),auxvec(1))
        call matmat2(auxmat(1,1),AECAderg(1,1,1),pizda(1,1))
        vv(1)=pizda(1,1)+pizda(2,2)
        vv(2)=pizda(2,1)-pizda(1,2)
@@@ -8490,12 -8177,12 +8510,12 @@@ C Cartesian derivatives
                s1=dip(4,jj,i)*dipderx(lll,kkk,4,kk,k)
              endif
  #endif
 -            call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
 +            call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
       &        auxvec(1))
 -            s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
 -            call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
 +            s2=0.5d0*scalar2(b1(1,k),auxvec(1))
 +            call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,l+1),
       &        auxvec(1))
 -            s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
 +            s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
              call matmat2(auxmat(1,1),AECAderx(1,1,lll,kkk,iii,1),
       &        pizda(1,1))
              vv(1)=pizda(1,1)+pizda(2,2)
@@@ -8534,7 -8221,7 +8554,7 @@@ c--------------------------------------
       & auxvec1(2),auxmat1(2,2)
        logical swap
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- C                                                                              C
+ C                                                                              C                       
  C      Parallel       Antiparallel                                             C
  C                                                                              C
  C          o             o                                                     C
@@@ -8542,10 -8229,10 +8562,10 @@@ C         /l\   /   \   /j
  C        /   \ /     \ /   \                                                   C
  C       /| o |o       o| o |\                                                  C
  C     \ j|/k\|      \  |/k\|l                                                  C
- C      \ /   \       \ /   \                                                   C
+ C      \ /   \       \ /   \                                                   C 
  C       o     \       o     \                                                  C
  C       i             i                                                        C
- C                                                                              C
+ C                                                                              C 
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  C
  C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
@@@ -8583,11 -8270,11 +8603,11 @@@ cd     & ' itl',itl,' itl1',itl
        call matvec2(AECA(1,1,imat),Ub2(1,k),auxvec(1))
        s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
        if (j.eq.l+1) then
 -        call matvec2(ADtEA1(1,1,3-imat),b1(1,itj1),auxvec1(1))
 -        s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
 +        call matvec2(ADtEA1(1,1,3-imat),b1(1,j+1),auxvec1(1))
 +        s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
        else
 -        call matvec2(ADtEA1(1,1,3-imat),b1(1,itl1),auxvec1(1))
 -        s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
 +        call matvec2(ADtEA1(1,1,3-imat),b1(1,l+1),auxvec1(1))
 +        s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
        endif
        call transpose2(EUg(1,1,k),auxmat(1,1))
        call matmat2(AECA(1,1,imat),auxmat(1,1),pizda(1,1))
@@@ -8611,11 -8298,11 +8631,11 @@@ C Derivatives in gamma(i-1
  #endif
          s2=0.5d0*scalar2(Ub2der(1,i),auxvec(1))
          if (j.eq.l+1) then
 -          call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itj1),auxvec1(1))
 -          s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
 +          call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,j+1),auxvec1(1))
 +          s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
          else
 -          call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itl1),auxvec1(1))
 -          s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
 +          call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,l+1),auxvec1(1))
 +          s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
          endif
          s4=0.25d0*scalar2(vv(1),Dtobr2der(1,i))
          if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then
@@@ -8644,11 -8331,11 +8664,11 @@@ C Derivatives in gamma(k-1
        call matvec2(AECA(1,1,imat),Ub2der(1,k),auxvec1(1))
        s2=0.5d0*scalar2(Ub2(1,i),auxvec1(1))
        if (j.eq.l+1) then
 -        call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itj1),auxvec1(1))
 -        s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
 +        call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,j+1),auxvec1(1))
 +        s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
        else
 -        call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itl1),auxvec1(1))
 -        s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
 +        call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,l+1),auxvec1(1))
 +        s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
        endif
        call transpose2(EUgder(1,1,k),auxmat1(1,1))
        call matmat2(AECA(1,1,imat),auxmat1(1,1),pizda(1,1))
@@@ -8714,12 -8401,12 +8734,12 @@@ C Cartesian derivatives
              s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
              if (j.eq.l+1) then
                call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
 -     &          b1(1,itj1),auxvec(1))
 -              s3=-0.5d0*scalar2(b1(1,itj),auxvec(1))
 +     &          b1(1,j+1),auxvec(1))
 +              s3=-0.5d0*scalar2(b1(1,j),auxvec(1))
              else
                call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
 -     &          b1(1,itl1),auxvec(1))
 -              s3=-0.5d0*scalar2(b1(1,itl),auxvec(1))
 +     &          b1(1,l+1),auxvec(1))
 +              s3=-0.5d0*scalar2(b1(1,l),auxvec(1))
              endif
              call matmat2(AECAderx(1,1,lll,kkk,iii,imat),auxmat(1,1),
       &        pizda(1,1))
@@@ -8819,12 -8506,12 +8839,12 @@@ cd      write (2,*) 'eello6_5',eello6_
  #ifdef MOMENT
        call transpose2(AEA(1,1,1),auxmat(1,1))
        call matmat2(EUg(1,1,i+1),auxmat(1,1),auxmat(1,1))
 -      ss1=scalar2(Ub2(1,i+2),b1(1,itl))
 +      ss1=scalar2(Ub2(1,i+2),b1(1,l))
        s1 = (auxmat(1,1)+auxmat(2,2))*ss1
  #endif
 -      call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
 +      call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1))
        call matvec2(AEA(1,1,1),vtemp1(1),vtemp1(1))
 -      s2 = scalar2(b1(1,itk),vtemp1(1))
 +      s2 = scalar2(b1(1,k),vtemp1(1))
  #ifdef MOMENT
        call transpose2(AEA(1,1,2),atemp(1,1))
        call matmat2(atemp(1,1),EUg(1,1,i+4),atemp(1,1))
@@@ -8839,7 -8526,7 +8859,7 @@@
        call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtemp(1,1))
        call matmat2(gtemp(1,1),EUg(1,1,i+3),gtemp(1,1)) 
        call matvec2(a_chuj(1,1,jj,i),Ub2(1,i+4),vtemp4(1)) 
 -      ss13 = scalar2(b1(1,itk),vtemp4(1))
 +      ss13 = scalar2(b1(1,k),vtemp4(1))
        s13 = (gtemp(1,1)+gtemp(2,2))*ss13
  #endif
  c      write (2,*) 's1,s2,s8,s12,s13',s1,s2,s8,s12,s13
@@@ -8873,12 -8560,12 +8893,12 @@@ C Derivatives in gamma(i+3
  #ifdef MOMENT
        call transpose2(AEA(1,1,1),auxmatd(1,1))
        call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
 -      ss1d=scalar2(Ub2der(1,i+2),b1(1,itl))
 +      ss1d=scalar2(Ub2der(1,i+2),b1(1,l))
        s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1d
  #endif
 -      call matvec2(EUgder(1,1,i+2),b1(1,itl),vtemp1d(1))
 +      call matvec2(EUgder(1,1,i+2),b1(1,l),vtemp1d(1))
        call matvec2(AEA(1,1,1),vtemp1d(1),vtemp1d(1))
 -      s2d = scalar2(b1(1,itk),vtemp1d(1))
 +      s2d = scalar2(b1(1,k),vtemp1d(1))
  #ifdef MOMENT
        call matvec2(Ug2der(1,1,i+2),dd(1,1,itk1),vtemp2d(1))
        s8d = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2d(1))
@@@ -8926,9 -8613,9 +8946,9 @@@ C Derivatives in gamma(i+5
        call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
        s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
  #endif
 -      call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1d(1))
 +      call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1d(1))
        call matvec2(AEAderg(1,1,1),vtemp1d(1),vtemp1d(1))
 -      s2d = scalar2(b1(1,itk),vtemp1d(1))
 +      s2d = scalar2(b1(1,k),vtemp1d(1))
  #ifdef MOMENT
        call transpose2(AEA(1,1,2),atempd(1,1))
        call matmat2(atempd(1,1),EUgder(1,1,i+4),atempd(1,1))
@@@ -8938,7 -8625,7 +8958,7 @@@
        s12d = scalar2(Ub2(1,i+2),vtemp3d(1))
  #ifdef MOMENT
        call matvec2(a_chuj(1,1,jj,i),Ub2der(1,i+4),vtemp4d(1)) 
 -      ss13d = scalar2(b1(1,itk),vtemp4d(1))
 +      ss13d = scalar2(b1(1,k),vtemp4d(1))
        s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
  #endif
  c      s1d=0.0d0
@@@ -8962,10 -8649,10 +8982,10 @@@ C Cartesian derivative
              call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
              s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
  #endif
 -            call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
 +            call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1))
              call matvec2(AEAderx(1,1,lll,kkk,iii,1),vtemp1(1),
       &          vtemp1d(1))
 -            s2d = scalar2(b1(1,itk),vtemp1d(1))
 +            s2d = scalar2(b1(1,k),vtemp1d(1))
  #ifdef MOMENT
              call transpose2(AEAderx(1,1,lll,kkk,iii,2),atempd(1,1))
              call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1))
@@@ -9009,7 -8696,7 +9029,7 @@@ c      s13d=0.0d
            derx_turn(lll,kkk,2) = derx_turn(lll,kkk,2)-0.5d0*s13d
            call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),Ub2(1,i+4),
       &      vtemp4d(1)) 
 -          ss13d = scalar2(b1(1,itk),vtemp4d(1))
 +          ss13d = scalar2(b1(1,k),vtemp4d(1))
            s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
            derx_turn(lll,kkk,1) = derx_turn(lll,kkk,1)-0.5d0*s13d
          enddo
diff --combined source/unres/src_MD-M/geomout.F
index c6a3944,6fd71c5..1e0c1c2
--- a/source/unres/src_MD-M/geomout.F
+++ b/source/unres/src_MD-M/geomout.F
@@@ -80,9 -80,15 +80,15 @@@ cmodel      write (iunit,'(a5,i6)') 'MO
  
        if (nss.gt.0) then
          do i=1,nss
+          if (dyn_ss) then
            write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') 
-      &         'SSBOND',i,'CYS',ihpb(i)-1-nres,
-      &                    'CYS',jhpb(i)-1-nres
+      &         'SSBOND',i,'CYS',idssb(i)-nnt+1,
+      &                    'CYS',jdssb(i)-nnt+1
+          else
+           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') 
+      &         'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,
+      &                    'CYS',jhpb(i)-nnt+1-nres
+          endif
          enddo
        endif
        
@@@ -91,7 -97,7 +97,7 @@@
        ires=0
        do i=nnt,nct
          iti=itype(i)
 -        if (iti.eq.21) then
 +        if (iti.eq.ntyp1) then
            ichain=ichain+1
            ires=0
            write (iunit,'(a)') 'TER'
@@@ -111,12 -117,12 +117,12 @@@
        enddo
        write (iunit,'(a)') 'TER'
        do i=nnt,nct-1
 -        if (itype(i).eq.21) cycle
 -        if (itype(i).eq.10 .and. itype(i+1).ne.21) then
 +        if (itype(i).eq.ntyp1) cycle
 +        if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
            write (iunit,30) ica(i),ica(i+1)
 -        else if (itype(i).ne.10 .and. itype(i+1).ne.21) then
 +        else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
            write (iunit,30) ica(i),ica(i+1),ica(i)+1
 -        else if (itype(i).ne.10 .and. itype(i+1).eq.21) then
 +        else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
            write (iunit,30) ica(i),ica(i)+1
          endif
        enddo
@@@ -124,7 -130,11 +130,11 @@@
          write (iunit,30) ica(nct),ica(nct)+1
        endif
        do i=1,nss
+        if (dyn_ss) then
+         write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
+        else
          write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
+        endif
        enddo
        write (iunit,'(a6)') 'ENDMDL'     
    10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
@@@ -200,7 -210,6 +210,7 @@@ c--------------------------------------
        include 'COMMON.INTERACT'
        include 'COMMON.NAMES'
        include 'COMMON.GEO'
 +      include 'COMMON.TORSION'
        write (iout,'(/a)') 'Geometry of the virtual chain.'
        write (iout,'(7a)') '  Res  ','         d','     Theta',
       & '       Phi','       Dsc','     Alpha','      Omega'
@@@ -279,8 -288,13 +289,13 @@@ c--------------------------------------
        open(icart,file=cartname,access="append")
  #endif
        write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
-       write (icart,'(i4,$)')
+       if (dyn_ss) then
+        write (icart,'(i4,$)')
+      &   nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)       
+       else
+        write (icart,'(i4,$)')
       &   nss,(ihpb(j),jhpb(j),j=1,nss)
+        endif
         write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,
       & (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),
       & (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
@@@ -322,8 -336,13 +337,13 @@@ c--------------------------------------
        call xdrffloat_(ixdrf, real(t_bath), iret)
        call xdrfint_(ixdrf, nss, iret) 
        do j=1,nss
+        if (dyn_ss) then
+         call xdrfint_(ixdrf, idssb(j)+nres, iret)
+         call xdrfint_(ixdrf, jdssb(j)+nres, iret)
+        else
          call xdrfint_(ixdrf, ihpb(j), iret)
          call xdrfint_(ixdrf, jhpb(j), iret)
+        endif
        enddo
        call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
        do i=1,nfrag
@@@ -346,8 -365,13 +366,13 @@@
        call xdrffloat(ixdrf, real(t_bath), iret)
        call xdrfint(ixdrf, nss, iret) 
        do j=1,nss
+        if (dyn_ss) then
+         call xdrfint(ixdrf, idssb(j)+nres, iret)
+         call xdrfint(ixdrf, jdssb(j)+nres, iret)
+        else
          call xdrfint(ixdrf, ihpb(j), iret)
          call xdrfint(ixdrf, jhpb(j), iret)
+        endif
        enddo
        call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
        do i=1,nfrag
diff --combined source/unres/src_MD-M/initialize_p.F
index a650241,26ba786..718e60e
--- a/source/unres/src_MD-M/initialize_p.F
+++ b/source/unres/src_MD-M/initialize_p.F
@@@ -80,9 -80,7 +80,9 @@@
        igeom=  8
        intin=  9
        ithep= 11
 +      ithep_pdb=51
        irotam=12
 +      irotam_pdb=52
        itorp= 13
        itordp= 23
        ielep= 14
@@@ -163,14 -161,10 +163,14 @@@ c      call memmon_print_usage(
  	rr0(i)=0.0D0
  	a0thet(i)=0.0D0
  	do j=1,2
 -	  athet(j,i)=0.0D0
 -	  bthet(j,i)=0.0D0
 +         do ichir1=-1,1
 +          do ichir2=-1,1
 +          athet(j,i,ichir1,ichir2)=0.0D0
 +          bthet(j,i,ichir1,ichir2)=0.0D0
 +          enddo
 +         enddo
          enddo
 -	do j=0,3
 +        do j=0,3
  	  polthet(j,i)=0.0D0
          enddo
  	do j=1,3
@@@ -194,39 -188,15 +194,39 @@@
        enddo
        nlob(ntyp1)=0
        dsc(ntyp1)=0.0D0
 -      do i=1,maxtor
 -	itortyp(i)=0
 -	do j=1,maxtor
 -	  do k=1,maxterm
 -	    v1(k,j,i)=0.0D0
 -	    v2(k,j,i)=0.0D0
 +      do i=-maxtor,maxtor
 +        itortyp(i)=0
 +cc      write (iout,*) "TU DOCHODZE",i,itortyp(i)
 +       do iblock=1,2
 +        do j=-maxtor,maxtor
 +          do k=1,maxterm
 +            v1(k,j,i,iblock)=0.0D0
 +            v2(k,j,i,iblock)=0.0D0
            enddo
          enddo
 +        enddo
        enddo
 +      do iblock=1,2
 +       do i=-maxtor,maxtor
 +        do j=-maxtor,maxtor
 +         do k=-maxtor,maxtor
 +          do l=1,maxtermd_1
 +            v1c(1,l,i,j,k,iblock)=0.0D0
 +            v1s(1,l,i,j,k,iblock)=0.0D0
 +            v1c(2,l,i,j,k,iblock)=0.0D0
 +            v1s(2,l,i,j,k,iblock)=0.0D0
 +          enddo !l
 +          do l=1,maxtermd_2
 +           do m=1,maxtermd_2
 +            v2c(m,l,i,j,k,iblock)=0.0D0
 +            v2s(m,l,i,j,k,iblock)=0.0D0
 +           enddo !m
 +          enddo !l
 +        enddo !k
 +       enddo !j
 +      enddo !i
 +      enddo !iblock
 +
        do i=1,maxres
  	itype(i)=0
  	itel(i)=0
@@@ -281,17 -251,11 +281,17 @@@ c--------------------------------------
        include 'COMMON.NAMES'
        include 'COMMON.FFIELD'
        data restyp /
 +     &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
 +     & 'DSG','DGN','DSN','DTH',
 +     &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
       &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 -     &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
 +     &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
 +     &'AIB','ABU','D'/
        data onelet /
 +     &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
 +     &'a','y','w','v','l','i','f','m','c','x',
       &'C','M','F','I','L','V','W','Y','A','G','T',
 -     &'S','Q','N','E','D','H','R','K','P','X'/
 +     &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
        data potname /'LJ','LJK','BP','GB','GBV'/
        data ename /
       &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
@@@ -377,6 -341,7 +377,7 @@@ cd    write (iout,*) 'ns=',ns,' nss=',n
  cd   &   (ihpb(i),jhpb(i),i=1,nss)
        do i=nnt,nct-1
          scheck=.false.
+         if (dyn_ss) goto 10
          do ii=1,nss
            if (ihpb(ii).eq.i+nres) then
              scheck=.true.
@@@ -596,9 -561,6 +597,9 @@@ C Partition local interaction
        iphi_end=iturn3_end+2
        iturn3_start=iturn3_start-1
        iturn3_end=iturn3_end-1
 +      call int_bounds(nres-3,itau_start,itau_end)
 +      itau_start=itau_start+3
 +      itau_end=itau_end+3
        call int_bounds(nres-3,iphi1_start,iphi1_end)
        iphi1_start=iphi1_start+3
        iphi1_end=iphi1_end+3
@@@ -1127,8 -1089,6 +1128,8 @@@ c        write (iout,*) "MPI_ROTAT2",MP
        idihconstr_end=ndih_constr
        iphid_start=iphi_start
        iphid_end=iphi_end-1
 +      itau_start=4
 +      itau_end=nres
        ibond_start=2
        ibond_end=nres-1
        ibondp_start=nnt
diff --combined source/unres/src_MD-M/readrtns_CSA.F
index 265b705,c2d0887..06a2146
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@@ -295,11 -295,11 +295,11 @@@ cd       endi
            do i=1,nrep
             iremd_m_total=iremd_m_total+remd_m(i)
            enddo
-           write (iout,*) 'Total number of replicas ',iremd_m_total
+            write (iout,*) 'Total number of replicas ',iremd_m_total
+           endif
           endif
-       endif
        if(me.eq.king.or..not.out1file) 
-      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
+      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
        return
        end
  c--------------------------------------------------------------------------
@@@ -345,6 -345,7 +345,7 @@@
        rest = index(controlcard,"REST").gt.0
        tbf = index(controlcard,"TBF").gt.0
        usampl = index(controlcard,"USAMPL").gt.0
+ 
        mdpdb = index(controlcard,"MDPDB").gt.0
        call reada(controlcard,"T_BATH",t_bath,300.0d0)
        call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1) 
@@@ -550,54 -551,54 +551,54 @@@
  C Body
  C
  C Read weights of the subsequent energy terms.
-       call card_concat(weightcard)
-       call reada(weightcard,'WLONG',wlong,1.0D0)
-       call reada(weightcard,'WSC',wsc,wlong)
-       call reada(weightcard,'WSCP',wscp,wlong)
-       call reada(weightcard,'WELEC',welec,1.0D0)
-       call reada(weightcard,'WVDWPP',wvdwpp,welec)
-       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
-       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
-       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
-       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
-       call reada(weightcard,'WTURN3',wturn3,1.0D0)
-       call reada(weightcard,'WTURN4',wturn4,1.0D0)
-       call reada(weightcard,'WTURN6',wturn6,1.0D0)
-       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
-       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
-       call reada(weightcard,'WBOND',wbond,1.0D0)
-       call reada(weightcard,'WTOR',wtor,1.0D0)
-       call reada(weightcard,'WTORD',wtor_d,1.0D0)
-       call reada(weightcard,'WANG',wang,1.0D0)
-       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
-       call reada(weightcard,'SCAL14',scal14,0.4D0)
-       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
-       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
-       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
-       call reada(weightcard,'TEMP0',temp0,300.0d0)
-       if (index(weightcard,'SOFT').gt.0) ipot=6
+        call card_concat(weightcard)
+        call reada(weightcard,'WLONG',wlong,1.0D0)
+        call reada(weightcard,'WSC',wsc,wlong)
+        call reada(weightcard,'WSCP',wscp,wlong)
+        call reada(weightcard,'WELEC',welec,1.0D0)
+        call reada(weightcard,'WVDWPP',wvdwpp,welec)
+        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
+        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
+        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
+        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
+        call reada(weightcard,'WTURN3',wturn3,1.0D0)
+        call reada(weightcard,'WTURN4',wturn4,1.0D0)
+        call reada(weightcard,'WTURN6',wturn6,1.0D0)
+        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
+        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
+        call reada(weightcard,'WBOND',wbond,1.0D0)
+        call reada(weightcard,'WTOR',wtor,1.0D0)
+        call reada(weightcard,'WTORD',wtor_d,1.0D0)
+        call reada(weightcard,'WANG',wang,1.0D0)
+        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+        call reada(weightcard,'SCAL14',scal14,0.4D0)
+        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
+        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
+        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
+        call reada(weightcard,'TEMP0',temp0,300.0d0)
+        if (index(weightcard,'SOFT').gt.0) ipot=6
  C 12/1/95 Added weight for the multi-body term WCORR
-       call reada(weightcard,'WCORRH',wcorr,1.0D0)
-       if (wcorr4.gt.0.0d0) wcorr=wcorr4
-       weights(1)=wsc
-       weights(2)=wscp
-       weights(3)=welec
-       weights(4)=wcorr
-       weights(5)=wcorr5
-       weights(6)=wcorr6
-       weights(7)=wel_loc
-       weights(8)=wturn3
-       weights(9)=wturn4
-       weights(10)=wturn6
-       weights(11)=wang
-       weights(12)=wscloc
-       weights(13)=wtor
-       weights(14)=wtor_d
-       weights(15)=wstrain
-       weights(16)=wvdwpp
-       weights(17)=wbond
-       weights(18)=scal14
-       weights(21)=wsccor
+        call reada(weightcard,'WCORRH',wcorr,1.0D0)
+        if (wcorr4.gt.0.0d0) wcorr=wcorr4
+        weights(1)=wsc
+        weights(2)=wscp
+        weights(3)=welec
+        weights(4)=wcorr
+        weights(5)=wcorr5
+        weights(6)=wcorr6
+        weights(7)=wel_loc
+        weights(8)=wturn3
+        weights(9)=wturn4
+        weights(10)=wturn6
+        weights(11)=wang
+        weights(12)=wscloc
+        weights(13)=wtor
+        weights(14)=wtor_d
+        weights(15)=wstrain
+        weights(16)=wvdwpp
+        weights(17)=wbond
+        weights(18)=scal14
+        weights(21)=wsccor
        if(me.eq.king.or..not.out1file)
       & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
       &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
@@@ -638,7 -639,7 +639,7 @@@
       &  'General scaling factor of SC-p interactions:',scalscp
        endif
        r0_corr=cutoff_corr-delt_corr
 -      do i=1,20
 +      do i=1,ntyp
          aad(i,1)=scalscp*aad(i,1)
          aad(i,2)=scalscp*aad(i,2)
          bad(i,1)=scalscp*bad(i,1)
@@@ -679,13 -680,37 +680,37 @@@
        call reada(weightcard,"V2SS",v2ss,7.61d0)
        call reada(weightcard,"V3SS",v3ss,13.7d0)
        call reada(weightcard,"EBR",ebr,-5.50D0)
+       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
+       do i=1,maxres
+         dyn_ss_mask(i)=.false.
+       enddo
+       do i=1,maxres-1
+         do j=i+1,maxres
+           dyn_ssbond_ij(i,j)=1.0d300
+         enddo
+       enddo
+       call reada(weightcard,"HT",Ht,0.0D0)
+       if (dyn_ss) then
+         ss_depth=ebr/wsc-0.25*eps(1,1)
+         Ht=Ht/wsc-0.25*eps(1,1)
+         akcm=akcm*wstrain/wsc
+         akth=akth*wstrain/wsc
+         akct=akct*wstrain/wsc
+         v1ss=v1ss*wstrain/wsc
+         v2ss=v2ss*wstrain/wsc
+         v3ss=v3ss*wstrain/wsc
+       else
+         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+       endif
+ 
        if(me.eq.king.or..not.out1file) then
         write (iout,*) "Parameters of the SS-bond potential:"
         write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
       & " AKCT",akct
         write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
-        write (iout,*) "EBR",ebr
- c       print *,'indpdb=',indpdb,' pdbref=',pdbref
+        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
+        write (iout,*)" HT",Ht
+        print *,'indpdb=',indpdb,' pdbref=',pdbref
        endif
        if (indpdb.gt.0 .or. pdbref) then
          read(inp,'(a)') pdbfile
@@@ -717,7 -742,7 +742,7 @@@ c        print *,'Finished reading pdb 
           maxsi=1000
           do i=2,nres-1
            iti=itype(i)
 -          if (iti.ne.10 .and. itype(i).ne.21) then
 +          if (iti.ne.10 .and. itype(i).ne.ntyp1) then
              nsi=0
              fail=.true.
              do while (fail.and.nsi.le.maxsi)
@@@ -749,8 -774,8 +774,8 @@@ C Assign initial virtual bond length
            vbld_inv(i)=vblinv
          enddo
          do i=2,nres-1
 -          vbld(i+nres)=dsc(itype(i))
 -          vbld_inv(i+nres)=dsc_inv(itype(i))
 +          vbld(i+nres)=dsc(iabs(itype(i)))
 +          vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
  c          write (iout,*) "i",i," itype",itype(i),
  c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
          enddo
@@@ -759,15 -784,15 +784,15 @@@ c      print *,nre
  c      print '(20i4)',(itype(i),i=1,nres)
        do i=1,nres
  #ifdef PROCOR
 -        if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 +        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
  #else
 -        if (itype(i).eq.21) then
 +        if (itype(i).eq.ntyp1) then
  #endif
            itel(i)=0
  #ifdef PROCOR
 -        else if (itype(i+1).ne.20) then
 +        else if (iabs(itype(i+1)).ne.20) then
  #else
 -        else if (itype(i).ne.20) then
 +        else if (iabs(itype(i)).ne.20) then
  #endif
  	  itel(i)=1
          else
@@@ -824,8 -849,8 +849,8 @@@ C 8/13/98 Set limits to generating the 
  #endif
        nct=nres
  cd      print *,'NNT=',NNT,' NCT=',NCT
 -      if (itype(1).eq.21) nnt=2
 -      if (itype(nres).eq.21) nct=nct-1
 +      if (itype(1).eq.ntyp1) nnt=2
 +      if (itype(nres).eq.ntyp1) nct=nct-1
        if (pdbref) then
          if(me.eq.king.or..not.out1file)
       &   write (iout,'(a,i3)') 'nsup=',nsup
@@@ -942,7 -967,7 +967,7 @@@ C initial geometry
                enddo
              enddo
              do i=nnt,nct
 -              if (itype(i).ne.10 .and. itype(i).ne.21) then
 +              if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
                  do j=1,3
                    dc(j,i+nres)=c(j,i+nres)-c(j,i) 
                    dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
@@@ -974,7 -999,6 +999,7 @@@
           enddo
           do i=2,nres-1
            omeg(i)=-120d0*deg2rad
 +          if (itype(i).le.0) omeg(i)=-omeg(i)
           enddo
          else
            if(me.eq.king.or..not.out1file)
@@@ -1012,18 -1036,7 +1037,7 @@@
     40       continue
            endif
  #else
-           do itrial=1,100
-             itmp=1
-             call gen_rand_conf(itmp,*30)
-             goto 40
-    30       write (iout,*) 'Failed to generate random conformation',
-      &        ', itrial=',itrial
-             write (*,*) 'Failed to generate random conformation',
-      &        ', itrial=',itrial
-           enddo
-           write (iout,'(a,i3,a)') 'Processor:',me,
-      &      ' error in generating random conformation.'
-           write (*,'(a,i3,a)') 'Processor:',me,
+           write (*,'(a)') 
       &      ' error in generating random conformation.'
            stop
     40     continue
@@@ -1055,18 -1068,35 +1069,35 @@@ C Generate distance constraints, if th
          write (iout,'(/a,i3,a)') 
       &  'The chain contains',ns,' disulfide-bridging cysteines.'
          write (iout,'(20i4)') (iss(i),i=1,ns)
+        if (dyn_ss) then
+           write(iout,*)"Running with dynamic disulfide-bond formation"
+        else
          write (iout,'(/a/)') 'Pre-formed links are:' 
  	do i=1,nss
  	  i1=ihpb(i)-nres
  	  i2=jhpb(i)-nres
  	  it1=itype(i1)
  	  it2=itype(i2)
- 	  if (me.eq.king.or..not.out1file)
-      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
+           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
       &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
       &    ebr,forcon(i)
  	enddo
  	write (iout,'(a)')
+        endif
+       endif
+       if (ns.gt.0.and.dyn_ss) then
+           do i=nss+1,nhpb
+             ihpb(i-nss)=ihpb(i)
+             jhpb(i-nss)=jhpb(i)
+             forcon(i-nss)=forcon(i)
+             dhpb(i-nss)=dhpb(i)
+           enddo
+           nhpb=nhpb-nss
+           nss=0
+           call hpb_partition
+           do i=1,ns
+             dyn_ss_mask(iss(i))=.true.
+           enddo
        endif
        if (i2ndstr.gt.0) call secstrp2dihc
  c      call geom_to_var(nvar,x)
@@@ -1123,17 -1153,19 +1154,19 @@@ C Read information about disulfide brid
        include 'COMMON.SETUP'
  C Read bridging residues.
        read (inp,*) ns,(iss(i),i=1,ns)
- c      print *,'ns=',ns
+       print *,'ns=',ns
        if(me.eq.king.or..not.out1file)
       &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
  C Check whether the specified bridging residues are cystines.
        do i=1,ns
  	if (itype(iss(i)).ne.1) then
  	  if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') 
-      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
+      &   'Do you REALLY think that the residue ',
+      &    restyp(itype(iss(i))),i,
       &   ' can form a disulfide bridge?!!!'
  	  write (*,'(2a,i3,a)') 
-      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
+      &   'Do you REALLY think that the residue ',
+      &    restyp(itype(iss(i))),i,
       &   ' can form a disulfide bridge?!!!'
  #ifdef MPI
  	 call MPI_Finalize(MPI_COMM_WORLD,ierror)
@@@ -1144,7 -1176,8 +1177,8 @@@
  C Read preformed bridges.
        if (ns.gt.0) then
        read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
-       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
+       if(fg_rank.eq.0)
+      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
        if (nss.gt.0) then
          nhpb=nss
  C Check if the residues involved in bridges are in the specified list of
@@@ -1212,7 -1245,7 +1246,7 @@@
          enddo
        enddo
        do i=nnt,nct
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
            do j=1,3
              dc(j,i+nres)=c(j,i+nres)-c(j,i)
              dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
@@@ -1292,7 -1325,7 +1326,7 @@@ C Set up variable list
        nvar=ntheta+nphi
        nside=0
        do i=2,nres-1
 -        if (itype(i).ne.10 .and. itype(i).ne.21) then
 +        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
  	  nside=nside+1
            ialph(i,1)=nvar+nside
  	  ialph(nside,2)=i
@@@ -1916,22 -1949,12 +1950,22 @@@ C Get parameter filenames and open the 
        open (itordp,file=tordname,status='old',readonly)
        call getenv_loc('SCCORPAR',sccorname)
        open (isccor,file=sccorname,status='old',readonly)
 +#ifndef CRYST_THETA
 +      call getenv_loc('THETPARPDB',thetname_pdb)
 +      print *,"thetname_pdb ",thetname_pdb
 +      open (ithep_pdb,file=thetname_pdb,status='old',action='read')
 +      print *,ithep_pdb," opened"
 +#endif
        call getenv_loc('FOURIER',fouriername)
        open (ifourier,file=fouriername,status='old',readonly)
        call getenv_loc('ELEPAR',elename)
        open (ielep,file=elename,status='old',readonly)
        call getenv_loc('SIDEPAR',sidename)
        open (isidep,file=sidename,status='old',readonly)
 +#ifndef CRYST_SC
 +      call getenv_loc('ROTPARPDB',rotname_pdb)
 +      open (irotam_pdb,file=rotname_pdb,status='old',action='read')
 +#endif
  #endif
  #ifndef OLDSCP
  C
@@@ -2002,7 -2025,7 +2036,7 @@@ c      print *,"Processor",myrank," fg_
        mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
        statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
        if (lentmp.gt.0)
-      &  call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
+      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
       &    //'.stat')
        rest2name=prefix(:ilen(prefix))//'.rst'
        if(usampl) then 
@@@ -2093,8 -2116,6 +2127,8 @@@ C Write file name
       &  thetname(:ilen(thetname))
        write (iout,*) "Rotamer parameter file          : ",
       &  rotname(:ilen(rotname))
 +      write (iout,*) "Thetpdb parameter file          : ",
 +     &  thetname_pdb(:ilen(thetname_pdb))
        write (iout,*) "Threading database              : ",
       &  patname(:ilen(patname))
        if (lentmp.ne.0) 
@@@ -2236,11 -2257,11 +2270,11 @@@ c        write (iout,*) i,ifrag_(1,i),i
  c            write (iout,*) "j",j," k",k
              ddjk=dist(j,k)
              if (constr_dist.eq.1) then
-               nhpb=nhpb+1
-               ihpb(nhpb)=j
-               jhpb(nhpb)=k
+             nhpb=nhpb+1
+             ihpb(nhpb)=j
+             jhpb(nhpb)=k
                dhpb(nhpb)=ddjk
-               forcon(nhpb)=wfrag_(i) 
+             forcon(nhpb)=wfrag_(i) 
              else if (constr_dist.eq.2) then
                if (ddjk.le.dist_cut) then
                  nhpb=nhpb+1
@@@ -2418,7 -2439,7 +2452,7 @@@
  #endif
        if (OKRandom) then
  c        r1 = prng_next(me)
-          r1=ran_number(0.0D0,1.0D0)
+         r1=ran_number(0.0D0,1.0D0)
          if(me.eq.king)
       &   write (iout,*) 'ran_num',r1
          if (r1.lt.0.0d0) OKRandom=.false.
diff --combined source/unres/src_MD-M/sc_move.F
index 2082e98,b6837fd..274767b
--- a/source/unres/src_MD-M/sc_move.F
+++ b/source/unres/src_MD-M/sc_move.F
@@@ -15,7 -15,9 +15,9 @@@ crc      implicit non
  c     Includes
        implicit real*8 (a-h,o-z)
        include 'DIMENSIONS'
+ #ifdef MPI
        include 'mpif.h'
+ #endif
        include 'COMMON.GEO'
        include 'COMMON.VAR'
        include 'COMMON.HEADER'
@@@ -211,7 -213,7 +213,7 @@@ c     Define what is meant by "neighbou
  
  c     Don't do glycine or ends
        i=itype(res_pick)
 -      if (i.eq.10 .or. i.eq.21) return
 +      if (i.eq.10 .or. i.eq.ntyp1) return
  
  c     Freeze everything (later will relax only selected side-chains)
        mask_r=.true.
@@@ -253,7 -255,7 +255,7 @@@ cd      print *,'new       ',(energy(k)
        n_try=0
        do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
  c     Move the selected residue (don't worry if it fails)
 -        call gen_side(itype(res_pick),theta(res_pick+1),
 +        call gen_side(iabs(itype(res_pick)),theta(res_pick+1),
       +       alph(res_pick),omeg(res_pick),fail)
  
  c     Minimize the side-chains starting from the new arrangement
@@@ -717,8 -719,8 +719,8 @@@ c     if (icall.eq.0) lprn=.true
        do i=iatsc_s,iatsc_e
  
  
 -        itypi=itype(i)
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -733,7 -735,7 +735,7 @@@
            do j=istart(i,iint),iend(i,iint)
            IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
              dscj_inv=dsc_inv(itypj)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
diff --combined source/unres/src_MD/CMakeLists.txt
index dea499c,ea7776c..8d21134
--- a/source/unres/src_MD/CMakeLists.txt
+++ b/source/unres/src_MD/CMakeLists.txt
@@@ -59,6 -59,7 +59,7 @@@ set(UNRES_MD_SRC
  	parmread.F 
  	pinorm.f 
  	printmat.f 
+ 	prng_32.F
  	q_measure.F 
  	randgens.f 
  	rattle.F 
@@@ -81,19 -82,6 +82,6 @@@
  	ssMD.F
  )
  
- if(Fortran_COMPILER_NAME STREQUAL "ifort")
-   set(UNRES_MD_SRC0 ${UNRES_MD_SRC0} prng.f ) 
- elseif(Fortran_COMPILER_NAME STREQUAL "mpif90")
-   set(UNRES_MD_SRC0 ${UNRES_MD_SRC0} prng.f )
- elseif(Fortran_COMPILER_NAME STREQUAL "f95")
-   set(UNRES_MD_SRC0 ${UNRES_MD_SRC0} prng.f )
- elseif(Fortran_COMPILER_NAME STREQUAL "gfortran")
-   set(UNRES_MD_SRC0 ${UNRES_MD_SRC0} prng.f )
- else()
-   set(UNRES_MD_SRC0 ${UNRES_MD_SRC0} prng_32.F )
- endif (Fortran_COMPILER_NAME STREQUAL "ifort")
- 
- 
  set(UNRES_MD_SRC3 
  	energy_p_new_barrier.F 
  	energy_p_new-sep_barrier.F 
@@@ -124,6 -112,7 +112,7 @@@ set(UNRES_MD_PP_SR
  	MP.F 
  	MREMD.F 
  	parmread.F 
+ 	prng_32.F
  	q_measure1.F 
  	q_measure3.F 
  	q_measure.F
@@@ -141,11 -130,6 +130,6 @@@
  	proc_proc.c 
  ) 
  
- 
- if(NOT Fortran_COMPILER_NAME STREQUAL "ifort")
-   set(UNRES_MD_PP_SRC ${UNRES_MD_PP_SRC} prng_32.F) 
- endif(NOT Fortran_COMPILER_NAME STREQUAL "ifort")
- 
  #================================================
  # Set comipiler flags for different sourcefiles  
  #================================================
@@@ -161,20 -145,15 +145,20 @@@ elseif (Fortran_COMPILER_NAME STREQUAL 
    set(FFLAGS2 "-std=legacy -I. ")
    #set(FFLAGS3 "-c -w -O3 -ipo -ipo_obj -opt_report" )
    set(FFLAGS3 "-std=legacy -I. " )
 +elseif (Fortran_COMPILER_NAME STREQUAL "g77")
 +  set(FFLAGS0 "-I ${CMAKE_CURRENT_SOURCE_DIR} " ) 
 +  set(FFLAGS1 "-g -I ${CMAKE_CURRENT_SOURCE_DIR} " ) 
 +  set(FFLAGS2 "-I ${CMAKE_CURRENT_SOURCE_DIR} ")
 +  set(FFLAGS3 "-I ${CMAKE_CURRENT_SOURCE_DIR} " )
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
  
  
  # Add MPI compiler flags
  if(UNRES_WITH_MPI)
-   set(FFLAGS0 "${FFLAGS0} -I${MPIF_INCLUDE_DIRECTORIES}")
-   set(FFLAGS1 "${FFLAGS1} -I${MPIF_INCLUDE_DIRECTORIES}")
-   set(FFLAGS2 "${FFLAGS2} -I${MPIF_INCLUDE_DIRECTORIES}")
-   set(FFLAGS3 "${FFLAGS3} -I${MPIF_INCLUDE_DIRECTORIES}")
+   set(FFLAGS0 "${FFLAGS0} -I${MPI_Fortran_INCLUDE_PATH}")
+   set(FFLAGS1 "${FFLAGS1} -I${MPI_Fortran_INCLUDE_PATH}")
+   set(FFLAGS2 "${FFLAGS2} -I${MPI_Fortran_INCLUDE_PATH}")
+   set(FFLAGS3 "${FFLAGS3} -I${MPI_Fortran_INCLUDE_PATH}")
  endif(UNRES_WITH_MPI)
  
  set_property(SOURCE ${UNRES_MD_SRC0} APPEND PROPERTY COMPILE_FLAGS ${FFLAGS0} )
@@@ -187,14 -166,14 +171,14 @@@ set_property(SOURCE ${UNRES_MD_SRC3} PR
  #=========================================
  if(UNRES_MD_FF STREQUAL "GAB" )
    # set preprocesor flags   
 -  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC" )
 +  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB" )
  
  #=========================================
  #  Settings for E0LL2Y force field
  #=========================================
  elseif(UNRES_MD_FF STREQUAL "E0LL2Y")
    # set preprocesor flags   
 -  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0" )
 +  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DSCCORPDB" )
  endif(UNRES_MD_FF STREQUAL "GAB")
  
  #=========================================
@@@ -217,12 -196,6 +201,12 @@@ elseif (Fortran_COMPILER_NAME STREQUAL 
  elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
    # Add old gfortran flags
    set(CPPFLAGS "${CPPFLAGS} -DG77") 
 +elseif (Fortran_COMPILER_NAME STREQUAL "g77")
 +  # Add old gfortran flags
 +  set(CPPFLAGS "${CPPFLAGS} -DG77")
 +else(Fortran_COMPILER_NAME STREQUAL "ifort")
 +  # Default preprocessor flag
 +  set(CPPFLAGS "${CPPFLAGS} -DPGI")
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
  
  #=========================================
@@@ -233,6 -206,13 +217,13 @@@ if (UNRES_WITH_MPI
  endif(UNRES_WITH_MPI)
  
  #=========================================
+ # Add 64-bit specific preprocessor flags
+ #=========================================
+ if (architektura STREQUAL "64")
+   set(CPPFLAGS "${CPPFLAGS} -DAMD64")
+ endif (architektura STREQUAL "64")
+ 
+ #=========================================
  # Apply preprocesor flags to *.F files
  #=========================================
  set_property(SOURCE ${UNRES_MD_PP_SRC} PROPERTY COMPILE_DEFINITIONS ${CPPFLAGS} )  
@@@ -243,10 -223,10 +234,10 @@@
  #========================================
  if(UNRES_WITH_MPI) 
    # binary with mpi
-   set(UNRES_BIN "unres_${Fortran_COMPILER_NAME}_MPICH_${UNRES_MD_FF}.exe")
+   set(UNRES_BIN "unresMD_${Fortran_COMPILER_NAME}_MPICH_${UNRES_MD_FF}.exe")
  else(UNRES_WITH_MPI)
    # binary without mpi
-   set(UNRES_BIN "unres_${Fortran_COMPILER_NAME}_single_${UNRES_MD_FF}.exe")
+   set(UNRES_BIN "unresMD_${Fortran_COMPILER_NAME}_single_${UNRES_MD_FF}.exe")
  endif(UNRES_WITH_MPI)  
  
  #=========================================
@@@ -292,7 -272,7 +283,7 @@@ set(UNRES_MD_SRCS ${UNRES_MD_SRC0} ${UN
  #=========================================
  add_executable(UNRES_BIN-MD ${UNRES_MD_SRCS} )
  set_target_properties(UNRES_BIN-MD PROPERTIES OUTPUT_NAME ${UNRES_BIN})
- #set_property(TARGET ${UNRES_BIN} PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin/unres/MD )
+ set_property(TARGET UNRES_BIN-MD PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  #add_dependencies (${UNRES_BIN} ${UNRES_XDRFLIB})
  
  
@@@ -301,13 -281,18 +292,18 @@@
  #=========================================
  # link MPI library (libmpich.a)  
  if(UNRES_WITH_MPI)
-   target_link_libraries( UNRES_BIN-MD ${MPIF_LIBRARIES} )
+   target_link_libraries( UNRES_BIN-MD ${MPI_Fortran_LIBRARIES} )
  endif(UNRES_WITH_MPI)
  # link libxdrf.a 
  #message("UNRES_XDRFLIB=${UNRES_XDRFLIB}")
  target_link_libraries( UNRES_BIN-MD xdrf )
  
  #=========================================
+ # Install Path
+ #=========================================
+ install(TARGETS UNRES_BIN-MD DESTINATION ${CMAKE_INSTALL_PREFIX}) 
+ 
+ #=========================================
  # TESTS 
  #=========================================
  
@@@ -337,7 -322,7 +333,7 @@@ FILE(WRITE ${CMAKE_CURRENT_BINARY_DIR}/
  export POT=GB
  export PREFIX=ala10
  #-----------------------------------------------------------------------------
- UNRES_BIN=./${UNRES_BIN}
+ UNRES_BIN=${CMAKE_BINARY_DIR}/bin/${UNRES_BIN}
  #-----------------------------------------------------------------------------
  DD=${CMAKE_SOURCE_DIR}/PARAM
  export BONDPAR=$DD/bond.parm
diff --combined source/unres/src_MD/COMMON.SCCOR
index 154de36,8de6d3c..7952bd1
--- a/source/unres/src_MD/COMMON.SCCOR
+++ b/source/unres/src_MD/COMMON.SCCOR
@@@ -2,17 -2,16 +2,17 @@@ cc Parameters of the SCCOR ter
        double precision v1sccor,v2sccor,vlor1sccor,
       &                 vlor2sccor,vlor3sccor,gloc_sc,
       &                 dcostau,dsintau,dtauangle,dcosomicron,
-      &                 domicron
+      &                 domicron,v0sccor
        integer nterm_sccor,isccortyp,nsccortyp,nlor_sccor
 -      common/sccor/v1sccor(maxterm_sccor,3,20,20),
 -     &    v2sccor(maxterm_sccor,3,20,20),
 +      common/sccor/v1sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp),
 +     &    v2sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp),
 +     &    v0sccor(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp),
 +     &    nterm_sccor(-ntyp:ntyp,-ntyp:ntyp),isccortyp(-ntyp:ntyp),
 +     &    nsccortyp,
 +     &    nlor_sccor(-ntyp:ntyp,-ntyp:ntyp),
       &    vlor1sccor(maxterm_sccor,20,20),
       &    vlor2sccor(maxterm_sccor,20,20),
       &    vlor3sccor(maxterm_sccor,20,20),gloc_sc(3,0:maxres2,10),
 -     &    v0sccor(ntyp,ntyp),
       &    dcostau(3,3,3,maxres2),dsintau(3,3,3,maxres2),
       &    dtauangle(3,3,3,maxres2),dcosomicron(3,3,3,maxres2),
 -     &    domicron(3,3,3,maxres2),
 -     &    nterm_sccor(ntyp,ntyp),isccortyp(ntyp),nsccortyp,
 -     &    nlor_sccor(ntyp,ntyp)
 +     &    domicron(3,3,3,maxres2)
diff --combined source/unres/src_MD/energy_p_new_barrier.F
index bf547f4,4f753e4..2e0ac37
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@@ -1449,7 -1449,7 +1449,7 @@@
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
  c            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              chi1=chi(itypi,itypj)
@@@ -4569,18 -4569,6 +4569,18 @@@ c     write (*,'(a,i2)') 'EBEND ICG=',i
  C Zero the energy function and its derivative at 0 or pi.
          call splinthet(theta(i),0.5d0*delta,ss,ssd)
          it=itype(i-1)
 +        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))
 +         endif
          if (i.gt.3) then
  #ifdef OSF
  	  phii=phi(i)
@@@ -4614,27 -4602,15 +4614,27 @@@ C dependent on the adjacent virtual-bon
  C In following comments this theta will be referred to as t_c.
          thet_pred_mean=0.0d0
          do k=1,2
 -          athetk=athet(k,it)
 -          bthetk=bthet(k,it)
 +            athetk=athet(k,it,ichir1,ichir2)
 +            bthetk=bthet(k,it,ichir1,ichir2)
 +          if (it.eq.10) then
 +             athetk=athet(k,itype1,ichir11,ichir12)
 +             bthetk=bthet(k,itype2,ichir21,ichir22)
 +          endif
            thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
          enddo
          dthett=thet_pred_mean*ssd
          thet_pred_mean=thet_pred_mean*ss+a0thet(it)
  C Derivatives of the "mean" values in gamma1 and gamma2.
 -        dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss
 -        dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss
 +        dthetg1=(-athet(1,it,ichir1,ichir2)*y(2)
 +     &+athet(2,it,ichir1,ichir2)*y(1))*ss
 +         dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2)
 +     &          +bthet(2,it,ichir1,ichir2)*z(1))*ss
 +         if (it.eq.10) then
 +      dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2)
 +     &+athet(2,itype1,ichir11,ichir12)*y(1))*ss
 +        dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2)
 +     &         +bthet(2,itype2,ichir21,ichir22)*z(1))*ss
 +         endif
          if (theta(i).gt.pi-delta) then
            call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
       &         E_tc0)
@@@ -4832,9 -4808,6 +4832,9 @@@
            enddo 
          endif
          if (i.lt.nres) then
 +
 +        if (iabs(itype(i+1)).eq.20) iblock=2
 +        if (iabs(itype(i+1)).ne.20) iblock=1
  #ifdef OSF
            phii1=phi(i+1)
            if (phii1.ne.phii1) phii1=150.0
@@@ -4855,7 -4828,7 +4855,7 @@@
              sinph2(k)=0.0d0
            enddo
          endif  
 -        ethetai=aa0thet(ityp1,ityp2,ityp3)
 +         ethetai=aa0thet(ityp1,ityp2,ityp3,iblock)
          do k=1,ndouble
            do l=1,k-1
              ccl=cosph1(l)*cosph2(k-l)
@@@ -4877,12 -4850,11 +4877,12 @@@
          enddo
          endif
          do k=1,ntheterm
 -          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
 -          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
 +          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k)
 +          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock)
       &      *coskt(k)
            if (lprn)
 -     &    write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
 +     &    write (iout,*) "k",k,
 +     &    "aathet",aathet(k,ityp1,ityp2,ityp3,iblock),
       &     " ethetai",ethetai
          enddo
          if (lprn) then
@@@ -4901,24 -4873,24 +4901,24 @@@
          endif
          do m=1,ntheterm2
            do k=1,nsingle
 -            aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
 -     &         +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
 -     &         +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
 -     &         +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
 +            aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
 +     &         +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)
 +     &         +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)
 +     &         +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)
              ethetai=ethetai+sinkt(m)*aux
              dethetai=dethetai+0.5d0*m*aux*coskt(m)
              dephii=dephii+k*sinkt(m)*(
 -     &          ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
 -     &          bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
 +     &          ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)-
 +     &          bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k))
              dephii1=dephii1+k*sinkt(m)*(
 -     &          eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
 -     &          ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
 +     &          eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)-
 +     &          ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k))
              if (lprn)
       &      write (iout,*) "m",m," k",k," bbthet",
 -     &         bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
 -     &         ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
 -     &         ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
 -     &         eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
 +     &         bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet",
 +     &         ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
 +     &         ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
 +     &         eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
            enddo
          enddo
          if (lprn)
@@@ -4926,33 -4898,28 +4926,33 @@@
          do m=1,ntheterm3
            do k=2,ndouble
              do l=1,k-1
 -              aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
 +       aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
 +     & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+
 +     & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
 +     & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)
 +
                ethetai=ethetai+sinkt(m)*aux
                dethetai=dethetai+0.5d0*m*coskt(m)*aux
                dephii=dephii+l*sinkt(m)*(
 -     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
 +     & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)-
 +     &  ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
 +     &  ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
 +     &  ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
 +
                dephii1=dephii1+(k-l)*sinkt(m)*(
 -     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
 +     &-ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
 +     & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
 +     & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)-
 +     & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
 +
                if (lprn) then
                write (iout,*) "m",m," k",k," l",l," ffthet",
 -     &            ffthet(l,k,m,ityp1,ityp2,ityp3),
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3),
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
 +     &            ffthet(l,k,m,ityp1,ityp2,ityp3,iblock),
 +     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet",
 +     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock),
 +     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock),
 +     &            " ethetai",ethetai
 +
                write (iout,*) cosph1ph2(l,k)*sinkt(m),
       &            cosph1ph2(k,l)*sinkt(m),
       &            sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
@@@ -4961,9 -4928,11 +4961,11 @@@
            enddo
          enddo
  10      continue
-         if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
-      &   i,theta(i)*rad2deg,phii*rad2deg,
+ c        lprn1=.true.
+         if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)') 
+      &  'ebe', i,theta(i)*rad2deg,phii*rad2deg,
       &   phii1*rad2deg,ethetai
+ c        lprn1=.false.
          etheta=etheta+ethetai
          if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
          if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
@@@ -5304,7 -5273,7 +5306,7 @@@
          cosfac=dsqrt(cosfac2)
          sinfac2=0.5d0/(1.0d0-costtab(i+1))
          sinfac=dsqrt(sinfac2)
 -        it=itype(i)
 +        it=iabs(itype(i))
          if (it.eq.10) goto 1
  c
  C  Compute the axes of tghe local cartesian coordinates system; store in
@@@ -5322,7 -5291,7 +5324,7 @@@ C     &   dc_norm(3,i+nres
            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)
 +          z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
          enddo     
  c       write (2,*) "i",i
  c       write (2,*) "x_prime",(x_prime(j),j=1,3)
@@@ -5354,7 -5323,7 +5356,7 @@@
  C Compute the energy of the ith side cbain
  C
  c        write (2,*) "xx",xx," yy",yy," zz",zz
 -        it=itype(i)
 +        it=iabs(itype(i))
          do j = 1,65
            x(j) = sc_parmin(j,it) 
          enddo
@@@ -5362,7 -5331,7 +5364,7 @@@
  Cc diagnostics - remove later
          xx1 = dcos(alph(2))
          yy1 = dsin(alph(2))*dcos(omeg(2))
 -        zz1 = -dsin(alph(2))*dsin(omeg(2))
 +        zz1 = -dsign(1.0, dfloat(itype(i)))*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
@@@ -5532,10 -5501,8 +5534,10 @@@ c     &   (dC_norm(j,i-1),j=1,3)," vbld
           dZZ_Ci1(k)=0.0d0
           dZZ_Ci(k)=0.0d0
           do j=1,3
 -           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
 -           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
 +           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
 +     &     *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
 +           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)
 +     &     *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
           enddo
            
           dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
@@@ -5821,8 -5788,6 +5823,8 @@@ c     lprn=.true
        etors=0.0D0
        do i=iphi_start,iphi_end
        etors_ii=0.0D0
 +c        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
 +c     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
          itori=itortyp(itype(i-2))
          itori1=itortyp(itype(i-1))
          phii=phi(i)
@@@ -5916,8 -5881,6 +5918,8 @@@ C Set lprn=.true. for debuggin
  c     lprn=.true.
        etors_d=0.0D0
        do i=iphid_start,iphid_end
 +c        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
 +c     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
          itori=itortyp(itype(i-2))
          itori1=itortyp(itype(i-1))
          itori2=itortyp(itype(i))
@@@ -5990,11 -5953,10 +5992,11 @@@ c        amino-acid residues
  C Set lprn=.true. for debugging
        lprn=.false.
  c      lprn=.true.
 -c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
 +c     write (iout,*) "EBACK_SC_COR",itau_start,itau_end
        esccor=0.0D0
        do i=itau_start,itau_end
          esccor_ii=0.0D0
 +        if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
          isccori=isccortyp(itype(i-2))
          isccori1=isccortyp(itype(i-1))
          phii=phi(i)
@@@ -6013,16 -5975,14 +6015,16 @@@ c   2 = Ca...Ca...Ca...S
  c   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.21).or.
 -     &      (itype(i-1).eq.21)))
 +     &      (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.21)))
 +     &     .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.21)))) cycle  
 -        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle
 -        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21))
 +     &      (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))
       & cycle
          do j=1,nterm_sccor(isccori,isccori1)
            v1ij=v1sccor(j,intertyp,isccori,isccori1)
@@@ -6037,9 -5997,9 +6039,9 @@@ c        write (iout,*) "WTF",intertyp,
  c     &gloc_sc(intertyp,i-3,icg)
          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,
 -     &  (v1sccor(j,intertyp,itori,itori1),j=1,6)
 -     & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
 +     &  restyp(itype(i-2)),i-2,restyp(itype(i-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
         enddo !intertyp
        enddo
diff --combined source/unres/src_MIN/CMakeLists.txt
index 3d930f7,ed0bb6e..5d8442f
--- a/source/unres/src_MIN/CMakeLists.txt
+++ b/source/unres/src_MIN/CMakeLists.txt
@@@ -33,7 -33,6 +33,7 @@@ set(UNRES_MIN_SRC
  	printmat.f 
  	randgens.f 
  	readrtns_min.F
 +	refsys.f
  	rescode.f 
  	refsys.f
  	rmdd.f 
@@@ -209,41 -208,17 +209,17 @@@ set(UNRES_MIN_SRCS ${UNRES_MIN_SRC0} ${
  #=========================================
  # Build the binary
  #=========================================
- add_executable(UNRES_BIN-MIN ${UNRES_MIN_SRCS} )
- set_target_properties(UNRES_BIN-MIN PROPERTIES OUTPUT_NAME ${UNRES_BIN})
+ add_executable(UNRES_MIN_BIN ${UNRES_MIN_SRCS} )
+ set_target_properties(UNRES_MIN_BIN PROPERTIES OUTPUT_NAME ${UNRES_BIN})
  
  if (Fortran_COMPILER_NAME STREQUAL "ifort")
-   target_link_libraries (UNRES_BIN-MIN ${CMAKE_THREAD_LIBS_INIT})
+   target_link_libraries (UNRES_MIN_BIN ${CMAKE_THREAD_LIBS_INIT})
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
+ set_property(TARGET UNRES_MIN_BIN PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  
- #set_property(TARGET ${UNRES_BIN} PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin/unres/MD )
  
  #=========================================
- # TESTS 
+ # Install Path
  #=========================================
+ install(TARGETS UNRES_MIN_BIN DESTINATION ${CMAKE_INSTALL_PREFIX})
  
- #-- Copy all the data files from the test directory into the source directory
- #SET(UNRES_TEST_FILES
- #	ala10.inp
- #    )
- 
- #FOREACH (UNRES_TEST_FILE ${UNRES_TEST_FILES})
- #      SET (unres_test_dest "${CMAKE_CURRENT_BINARY_DIR}/${UNRES_TEST_FILE}")
- #      MESSAGE (STATUS " Copying ${UNRES_TEST_FILE} from ${CMAKE_SOURCE_DIR}/examples/unres/MD/ff_gab/${UNRES_TEST_FILE} to ${unres_test_dest}")
- #      ADD_CUSTOM_COMMAND (
- #          TARGET     ${UNRES_BIN}
- #          POST_BUILD
- #          COMMAND    ${CMAKE_COMMAND} -E copy ${CMAKE_SOURCE_DIR}/examples/unres/MD/ff_gab/${UNRES_TEST_FILE} ${unres_test_dest}
- #      )
- #ENDFOREACH (UNRES_TEST_FILE ${UNRES_TEST_FILES})
- 
- #=========================================
- # Generate data test files
- #=========================================
- 
- #if(NOT UNRES_WITH_MPI)
- 
- #  add_test(NAME UNRES_MD_Ala10 COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_single_ala.sh )
- 
- #endif(NOT UNRES_WITH_MPI)
-  
diff --combined source/wham/src-M/CMakeLists.txt
index 2576a15,afa0f19..312fb87
--- a/source/wham/src-M/CMakeLists.txt
+++ b/source/wham/src-M/CMakeLists.txt
@@@ -96,9 -96,11 +96,11 @@@ set(UNRES_WHAM_M_PP_SR
  # Set comipiler flags for different sourcefiles  
  #================================================
  if (Fortran_COMPILER_NAME STREQUAL "ifort")
-   set(FFLAGS0 "-g -CB -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres -I${MPIF_INCLUDE_DIRECTORIES}" ) 
+   set(FFLAGS0 "-g -CB -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" ) 
  elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
-   set(FFLAGS0 "-std=legacy -g -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres -I${MPIF_INCLUDE_DIRECTORIES}" ) 
+   set(FFLAGS0 "-std=legacy -g -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" ) 
+ else ()
+   set(FFLAGS0 "-g -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" ) 
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
  
  
@@@ -106,7 -108,7 +108,7 @@@
  # Add MPI compiler flags
  #=========================================
  if(UNRES_WITH_MPI)
-   set(FFLAGS0 "${FFLAGS0} -I${MPIF_INCLUDE_DIRECTORIES}")
+   set(FFLAGS0 "${FFLAGS0} -I${MPI_Fortran_INCLUDE_PATH}")
  endif(UNRES_WITH_MPI)
  
  set_property(SOURCE ${UNRES_WHAM_M_SRC0} PROPERTY COMPILE_FLAGS ${FFLAGS0} )
@@@ -118,7 -120,6 +120,7 @@@ if(UNRES_MD_FF STREQUAL "GAB" 
    # set preprocesor flags   
    set(CPPFLAGS "PROCOR  -DSPLITELE -DCRYST_BOND -DCRYST_THETA -DCRYST_SC  -DSCCORPDB" )
  
 +
  #=========================================
  #  Settings for E0LL2Y force field
  #=========================================
@@@ -152,9 -153,6 +154,9 @@@ elseif (Fortran_COMPILER_NAME STREQUAL 
  elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
    # Add old gfortran flags
    set(CPPFLAGS "${CPPFLAGS} -DG77") 
 +else (Fortran_COMPILER_NAME STREQUAL "ifort")
 +  # Default preprocessor flags
 +  set(CPPFLAGS "${CPPFLAGS} -DPGI")
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
  
  #=========================================
@@@ -178,7 -176,7 +180,7 @@@ set_property(SOURCE ${UNRES_WHAM_M_PP_S
  #========================================
  #  Setting binary name
  #========================================
- set(UNRES_WHAM_M_BIN "wham_${Fortran_COMPILER_NAME}.exe")
+ set(UNRES_WHAM_M_BIN "wham_M_${Fortran_COMPILER_NAME}.exe")
  
  #=========================================
  # cinfo.f workaround for CMake
@@@ -222,18 -220,24 +224,24 @@@ set(UNRES_WHAM_M_SRCS ${UNRES_WHAM_M_SR
  #=========================================
  add_executable(UNRES_WHAM_M_BIN ${UNRES_WHAM_M_SRCS} )
  set_target_properties(UNRES_WHAM_M_BIN PROPERTIES OUTPUT_NAME ${UNRES_WHAM_M_BIN})
- 
- #set_property(TARGET ${UNRES_BIN} PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin/unres/MD )
+ set_property(TARGET UNRES_WHAM_M_BIN PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  #add_dependencies (${UNRES_BIN} ${UNRES_XDRFLIB})
  
  #=========================================
  # Link libraries
  #=========================================
  # link MPI library (libmpich.a)  
- target_link_libraries( UNRES_WHAM_M_BIN ${MPIF_LIBRARIES} )
+ target_link_libraries( UNRES_WHAM_M_BIN ${MPI_Fortran_LIBRARIES} )
  # link libxdrf.a 
  target_link_libraries( UNRES_WHAM_M_BIN xdrf )
  
+ 
+ #=========================================
+ # Install Path
+ #=========================================
+ install(TARGETS UNRES_WHAM_M_BIN DESTINATION ${CMAKE_INSTALL_PREFIX})
+ 
+ 
  #=========================================
  # TESTS 
  #=========================================
diff --combined source/wham/src/CMakeLists.txt
index 6f762b1,e7f990e..43a8912
--- a/source/wham/src/CMakeLists.txt
+++ b/source/wham/src/CMakeLists.txt
@@@ -97,6 -97,8 +97,8 @@@ if (Fortran_COMPILER_NAME STREQUAL "ifo
    set(FFLAGS0 "-mcmodel=medium -g -CB -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" ) 
  elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
    set(FFLAGS0 "-std=legacy -g -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" ) 
+ else ()
+   set(FFLAGS0 "-g -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" )   
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
  
  
@@@ -104,7 -106,7 +106,7 @@@
  # Add MPI compiler flags
  #=========================================
  if(UNRES_WITH_MPI)
-   set(FFLAGS0 "${FFLAGS0} -I${MPIF_INCLUDE_DIRECTORIES}")
+   set(FFLAGS0 "${FFLAGS0} -I${MPI_Fortran_INCLUDE_PATH}")
  endif(UNRES_WITH_MPI)
  
  set_property(SOURCE ${UNRES_WHAM_SRC0} PROPERTY COMPILE_FLAGS ${FFLAGS0} )
@@@ -152,9 -154,6 +154,9 @@@ elseif (Fortran_COMPILER_NAME STREQUAL 
  elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
    # Add old gfortran flags
    set(CPPFLAGS "${CPPFLAGS} -DG77") 
 +else (Fortran_COMPILER_NAME STREQUAL "ifort")
 +  # Default preprocessor flags 
 +  set(CPPFLAGS "${CPPFLAGS} -DPGI")
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
  
  #=========================================
@@@ -222,19 -221,24 +224,24 @@@ set(UNRES_WHAM_SRCS ${UNRES_WHAM_SRC0} 
  #=========================================
  add_executable(UNRES_WHAM_BIN ${UNRES_WHAM_SRCS} )
  set_target_properties(UNRES_WHAM_BIN PROPERTIES OUTPUT_NAME ${UNRES_WHAM_BIN})
- 
- #set_property(TARGET ${UNRES_BIN} PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin/unres/MD )
+ set_property(TARGET UNRES_WHAM_BIN PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  #add_dependencies (${UNRES_BIN} ${UNRES_XDRFLIB})
  
  #=========================================
  # Link libraries
  #=========================================
- # link MPI library (libmpich.a)  
- target_link_libraries( UNRES_WHAM_BIN ${MPIF_LIBRARIES} )
+ # link MPI libraries
+ target_link_libraries( UNRES_WHAM_BIN ${MPI_Fortran_LIBRARIES} )
  # link libxdrf.a 
  target_link_libraries( UNRES_WHAM_BIN xdrf )
  
  #=========================================
+ # Install Path
+ #=========================================
+ install(TARGETS UNRES_WHAM_BIN DESTINATION ${CMAKE_INSTALL_PREFIX})
+ 
+ 
+ #=========================================
  # TESTS 
  #=========================================
  
diff --combined source/wham/src/energy_p_new.F
index a53e6e9,13fe796..7244671
--- a/source/wham/src/energy_p_new.F
+++ b/source/wham/src/energy_p_new.F
@@@ -229,7 -229,7 +229,7 @@@
       &   +wturn3*fact(2)*gel_loc_turn3(i)
       &   +wturn6*fact(5)*gel_loc_turn6(i)
       &   +wel_loc*fact(2)*gel_loc_loc(i)
 -     &   +wsccor*fact(1)*gsccor_loc(i)
 +c     &   +wsccor*fact(1)*gsccor_loc(i)
        enddo
        endif
        return
@@@ -369,9 -369,8 +369,9 @@@ cd    print *,'Entering ELJ nnt=',nnt,
        evdw=0.0D0
        evdw_t=0.0d0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        if (itypi.eq.ntyp1) cycle
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -384,8 -383,7 +384,8 @@@
  cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
  cd   &                  'iend=',iend(i,iint)
            do j=istart(i,iint),iend(i,iint)
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
 +            if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              zj=c(3,nres+j)-zi
@@@ -542,9 -540,8 +542,9 @@@ c     print *,'Entering ELJK nnt=',nnt,
        evdw=0.0D0
        evdw_t=0.0d0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        if (itypi.eq.ntyp1) cycle
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -553,8 -550,7 +553,8 @@@ C Calculate SC interaction energy
  C
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
 +            if (itypj.eq.ntyp1) cycle
              xj=c(1,nres+j)-xi
              yj=c(2,nres+j)-yi
              zj=c(3,nres+j)-zi
@@@ -655,8 -651,8 +655,8 @@@ c     els
  c     endif
        ind=0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -670,7 -666,7 +670,7 @@@
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
              dscj_inv=vbld_inv(j+nres)
              chi1=chi(itypi,itypj)
              chi2=chi(itypj,itypi)
@@@ -792,8 -788,8 +792,8 @@@ c     print *,'Entering EGB nnt=',nnt,
  c      if (icall.gt.0) lprn=.true.
        ind=0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -822,7 -818,7 +822,7 @@@ c              if (energy_dec) write (i
  c     &                        'evdw',i,j,evdwij,' ss'
              ELSE
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
              dscj_inv=vbld_inv(j+nres)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
@@@ -954,8 -950,8 +954,8 @@@ c     print *,'Entering EGB nnt=',nnt,
  c      if (icall.gt.0) lprn=.true.
        ind=0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        itypi1=itype(i+1)
 +        itypi=iabs(itype(i))
 +        itypi1=iabs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -969,7 -965,7 +969,7 @@@
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
              dscj_inv=vbld_inv(j+nres)
              sig0ij=sigma(itypi,itypj)
              r0ij=r0(itypi,itypj)
@@@ -2808,7 -2804,7 +2808,7 @@@ c     &   " iscp",(iscpstart(i,j),iscpe
          do iint=1,nscp_gr(i)
  
          do j=iscpstart(i,iint),iscpend(i,iint)
 -          itypj=itype(j)
 +          itypj=iabs(itype(j))
  C Uncomment following three lines for SC-p interactions
  c         xj=c(1,nres+j)-xi
  c         yj=c(2,nres+j)-yi
@@@ -2922,8 -2918,7 +2922,8 @@@ C 24/11/03 AL: SS bridges handled separ
  C    distance and angle dependent SS bond potential.
          if (.not.dyn_ss .and. i.le.nss) then
  C 15/02/13 CC dynamic SSbond - additional check
 -        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
 +         if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
 +     & iabs(itype(jjj)).eq.1) then
            call ssbond_ene(iii,jjj,eij)
            ehpb=ehpb+2*eij
           endif
@@@ -3024,7 -3019,7 +3024,7 @@@
        include 'COMMON.VAR'
        include 'COMMON.IOUNITS'
        double precision erij(3),dcosom1(3),dcosom2(3),gg(3)
 -      itypi=itype(i)
 +      itypi=iabs(itype(i))
        xi=c(1,nres+i)
        yi=c(2,nres+i)
        zi=c(3,nres+i)
@@@ -3032,7 -3027,7 +3032,7 @@@
        dyi=dc_norm(2,nres+i)
        dzi=dc_norm(3,nres+i)
        dsci_inv=dsc_inv(itypi)
 -      itypj=itype(j)
 +      itypj=iabs(itype(j))
        dscj_inv=dsc_inv(itypj)
        xj=c(1,nres+j)-xi
        yj=c(2,nres+j)-yi
@@@ -3124,7 -3119,7 +3124,7 @@@
  c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
  c
        do i=nnt,nct
 -        iti=itype(i)
 +        iti=iabs(itype(i))
          if (iti.ne.10) then
            nbi=nbondterm(iti)
            if (nbi.eq.1) then
@@@ -3206,18 -3201,6 +3206,18 @@@ c      write (iout,*) ithet_start,ithet
  C Zero the energy function and its derivative at 0 or pi.
          call splinthet(theta(i),0.5d0*delta,ss,ssd)
          it=itype(i-1)
 +        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))
 +         endif
  c        if (i.gt.ithet_start .and. 
  c     &     (itel(i-1).eq.0 .or. itel(i-2).eq.0)) goto 1215
  c        if (i.gt.3 .and. (i.le.4 .or. itel(i-3).ne.0)) then
@@@ -3273,12 -3256,8 +3273,12 @@@ C dependent on the adjacent virtual-bon
  C In following comments this theta will be referred to as t_c.
          thet_pred_mean=0.0d0
          do k=1,2
 -          athetk=athet(k,it)
 -          bthetk=bthet(k,it)
 +            athetk=athet(k,it,ichir1,ichir2)
 +            bthetk=bthet(k,it,ichir1,ichir2)
 +          if (it.eq.10) then
 +             athetk=athet(k,itype1,ichir11,ichir12)
 +             bthetk=bthet(k,itype2,ichir21,ichir22)
 +          endif
            thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
          enddo
  c        write (iout,*) "thet_pred_mean",thet_pred_mean
@@@ -3286,16 -3265,8 +3286,16 @@@
          thet_pred_mean=thet_pred_mean*ss+a0thet(it)
  c        write (iout,*) "thet_pred_mean",thet_pred_mean
  C Derivatives of the "mean" values in gamma1 and gamma2.
 -        dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss
 -        dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss
 +        dthetg1=(-athet(1,it,ichir1,ichir2)*y(2)
 +     &+athet(2,it,ichir1,ichir2)*y(1))*ss
 +         dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2)
 +     &          +bthet(2,it,ichir1,ichir2)*z(1))*ss
 +         if (it.eq.10) then
 +      dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2)
 +     &+athet(2,itype1,ichir11,ichir12)*y(1))*ss
 +        dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2)
 +     &         +bthet(2,itype2,ichir21,ichir22)*z(1))*ss
 +         endif
          if (theta(i).gt.pi-delta) then
            call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
       &         E_tc0)
@@@ -3466,13 -3437,11 +3466,13 @@@
        etheta=0.0D0
  c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
        do i=ithet_start,ithet_end
 +        if (iabs(itype(i+1)).eq.20) iblock=2
 +        if (iabs(itype(i+1)).ne.20) iblock=1
          dethetai=0.0d0
          dephii=0.0d0
          dephii1=0.0d0
          theti2=0.5d0*theta(i)
 -        ityp2=ithetyp(itype(i-1))
 +        ityp2=ithetyp((itype(i-1)))
          do k=1,nntheterm
            coskt(k)=dcos(k*theti2)
            sinkt(k)=dsin(k*theti2)
@@@ -3484,7 -3453,7 +3484,7 @@@
  #else
            phii=phi(i)
  #endif
 -          ityp1=ithetyp(itype(i-2))
 +          ityp1=ithetyp(iabs(itype(i-2)))
            do k=1,nsingle
              cosph1(k)=dcos(k*phii)
              sinph1(k)=dsin(k*phii)
@@@ -3505,7 -3474,7 +3505,7 @@@
  #else
            phii1=phi(i+1)
  #endif
 -          ityp3=ithetyp(itype(i))
 +          ityp3=ithetyp((itype(i)))
            do k=1,nsingle
              cosph2(k)=dcos(k*phii1)
              sinph2(k)=dsin(k*phii1)
@@@ -3521,7 -3490,7 +3521,7 @@@
  c        write (iout,*) "i",i," ityp1",itype(i-2),ityp1,
  c     &   " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3
  c        call flush(iout)
 -        ethetai=aa0thet(ityp1,ityp2,ityp3)
 +        ethetai=aa0thet(ityp1,ityp2,ityp3,iblock)
          do k=1,ndouble
            do l=1,k-1
              ccl=cosph1(l)*cosph2(k-l)
@@@ -3543,12 -3512,11 +3543,12 @@@
          enddo
          endif
          do k=1,ntheterm
 -          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
 -          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
 +          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k)
 +          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock)
       &      *coskt(k)
            if (lprn)
 -     &    write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
 +     &    write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3,
 +     &      iblock),
       &     " ethetai",ethetai
          enddo
          if (lprn) then
@@@ -3567,24 -3535,24 +3567,24 @@@
          endif
          do m=1,ntheterm2
            do k=1,nsingle
 -            aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
 -     &         +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
 -     &         +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
 -     &         +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
 +            aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
 +     &         +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)
 +     &         +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)
 +     &         +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)
              ethetai=ethetai+sinkt(m)*aux
              dethetai=dethetai+0.5d0*m*aux*coskt(m)
              dephii=dephii+k*sinkt(m)*(
 -     &          ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
 -     &          bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
 +     &          ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)-
 +     &          bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k))
              dephii1=dephii1+k*sinkt(m)*(
 -     &          eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
 -     &          ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
 +     &          eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)-
 +     &          ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k))
              if (lprn)
       &      write (iout,*) "m",m," k",k," bbthet",
 -     &         bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
 -     &         ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
 -     &         ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
 -     &         eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
 +     &         bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet",
 +     &         ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
 +     &         ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
 +     &         eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
            enddo
          enddo
          if (lprn)
@@@ -3592,29 -3560,28 +3592,29 @@@
          do m=1,ntheterm3
            do k=2,ndouble
              do l=1,k-1
 -              aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
 +              aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
 +     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+
 +     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
 +     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)
                ethetai=ethetai+sinkt(m)*aux
                dethetai=dethetai+0.5d0*m*coskt(m)*aux
                dephii=dephii+l*sinkt(m)*(
 -     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
 +     &           -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)-
 +     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
 +     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
 +     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
                dephii1=dephii1+(k-l)*sinkt(m)*(
 -     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
 +     &           -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
 +     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
 +     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)-
 +     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
                if (lprn) then
                write (iout,*) "m",m," k",k," l",l," ffthet",
 -     &            ffthet(l,k,m,ityp1,ityp2,ityp3),
 -     &            ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
 -     &            ggthet(l,k,m,ityp1,ityp2,ityp3),
 -     &            ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
 +     &            ffthet(l,k,m,ityp1,ityp2,ityp3,iblock),
 +     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet",
 +     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock),
 +     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ethetai",
 +     &            ethetai
                write (iout,*) cosph1ph2(l,k)*sinkt(m),
       &            cosph1ph2(k,l)*sinkt(m),
       &            sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
@@@ -3623,10 -3590,13 +3623,13 @@@
            enddo
          enddo
  10      continue
-         if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
-      &   i,theta(i)*rad2deg,phii*rad2deg,
+ c        lprn1=.true.
+         if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)') 
+      &  'ebe',i,theta(i)*rad2deg,phii*rad2deg,
       &   phii1*rad2deg,ethetai
+ c        lprn1=.false.
          etheta=etheta+ethetai
+         
          if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
          if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
          gloc(nphi+i-2,icg)=wang*dethetai
@@@ -3661,7 -3631,7 +3664,7 @@@ c     write (iout,'(a)') 'ESC
        do i=loc_start,loc_end
          it=itype(i)
          if (it.eq.10) goto 1
 -        nlobit=nlob(it)
 +        nlobit=nlob(iabs(it))
  c       print *,'i=',i,' it=',it,' nlobit=',nlobit
  c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
          theti=theta(i+1)-pipol
@@@ -3816,7 -3786,7 +3819,7 @@@ C Compute the contribution to SC energ
          do iii=-1,1
  
            do j=1,nlobit
 -            expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin)
 +            expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin)
  cd          print *,'j=',j,' expfac=',expfac
              escloc_i=escloc_i+expfac
              do k=1,3
@@@ -3897,7 -3867,7 +3900,7 @@@ C Compute the contribution to SC energ
  
        dersc12=0.0d0
        do j=1,nlobit
 -        expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin)
 +        expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin)
          escloc_i=escloc_i+expfac
          do k=1,2
            dersc(k)=dersc(k)+Ax(k,j)*expfac
@@@ -3960,7 -3930,7 +3963,7 @@@
          cosfac=dsqrt(cosfac2)
          sinfac2=0.5d0/(1.0d0-costtab(i+1))
          sinfac=dsqrt(sinfac2)
 -        it=itype(i)
 +        it=iabs(itype(i))
          if (it.eq.10) goto 1
  c
  C  Compute the axes of tghe local cartesian coordinates system; store in
@@@ -3978,7 -3948,7 +3981,7 @@@ C     &   dc_norm(3,i+nres
            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)
 +          z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
          enddo     
  c       write (2,*) "i",i
  c       write (2,*) "x_prime",(x_prime(j),j=1,3)
@@@ -4010,7 -3980,7 +4013,7 @@@
  C Compute the energy of the ith side cbain
  C
  c        write (2,*) "xx",xx," yy",yy," zz",zz
 -        it=itype(i)
 +        it=iabs(itype(i))
          do j = 1,65
            x(j) = sc_parmin(j,it) 
          enddo
@@@ -4018,7 -3988,7 +4021,7 @@@
  Cc diagnostics - remove later
          xx1 = dcos(alph(2))
          yy1 = dsin(alph(2))*dcos(omeg(2))
 -        zz1 = -dsin(alph(2))*dsin(omeg(2))
 +        zz1 = -dsign(1.0d0,itype(i))*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
@@@ -4190,11 -4160,8 +4193,11 @@@ c     &   (dC_norm(j,i-1),j=1,3)," vbld
           dZZ_Ci1(k)=0.0d0
           dZZ_Ci(k)=0.0d0
           do j=1,3
 -           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
 -           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
 +            dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
 +     & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
 +            dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)
 +     &  *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
 +
           enddo
            
           dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
@@@ -4430,19 -4397,14 +4433,19 @@@ c      lprn=.true
        etors=0.0D0
        do i=iphi_start,iphi_end
          if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215
 +         if (iabs(itype(i)).eq.20) then
 +         iblock=2
 +         else
 +         iblock=1
 +         endif
          itori=itortyp(itype(i-2))
          itori1=itortyp(itype(i-1))
          phii=phi(i)
          gloci=0.0D0
  C Regular cosine and sine terms
 -        do j=1,nterm(itori,itori1)
 -          v1ij=v1(j,itori,itori1)
 -          v2ij=v2(j,itori,itori1)
 +        do j=1,nterm(itori,itori1,iblock)
 +          v1ij=v1(j,itori,itori1,iblock)
 +          v2ij=v2(j,itori,itori1,iblock)
            cosphi=dcos(j*phii)
            sinphi=dsin(j*phii)
            etors=etors+v1ij*cosphi+v2ij*sinphi
@@@ -4455,7 -4417,7 +4458,7 @@@ C          [v2 cos(phi/2)+v3 sin(phi/2)
  C
          cosphi=dcos(0.5d0*phii)
          sinphi=dsin(0.5d0*phii)
 -        do j=1,nlor(itori,itori1)
 +        do j=1,nlor(itori,itori1,iblock)
            vl1ij=vlor1(j,itori,itori1)
            vl2ij=vlor2(j,itori,itori1)
            vl3ij=vlor3(j,itori,itori1)
@@@ -4466,11 -4428,11 +4469,11 @@@
            gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
          enddo
  C Subtract the constant term
 -        etors=etors-v0(itori,itori1)
 +        etors=etors-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,
 -     &  (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
 +     &  (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),j=1,6)
          gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci
  c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
   1215   continue
@@@ -4535,20 -4497,12 +4538,20 @@@ c     lprn=.true
          phii1=phi(i+1)
          gloci1=0.0D0
          gloci2=0.0D0
 +        iblock=1
 +        if (iabs(itype(i+1)).eq.20) iblock=2
  C Regular cosine and sine terms
 -        do j=1,ntermd_1(itori,itori1,itori2)
 -          v1cij=v1c(1,j,itori,itori1,itori2)
 -          v1sij=v1s(1,j,itori,itori1,itori2)
 -          v2cij=v1c(2,j,itori,itori1,itori2)
 -          v2sij=v1s(2,j,itori,itori1,itori2)
 +c c       do j=1,ntermd_1(itori,itori1,itori2,iblock)
 +c          v1cij=v1c(1,j,itori,itori1,itori2,iblock)
 +c          v1sij=v1s(1,j,itori,itori1,itori2,iblock)
 +c          v2cij=v1c(2,j,itori,itori1,itori2,iblock)
 +c          v2sij=v1s(2,j,itori,itori1,itori2,iblock)
 +       do j=1,ntermd_1(itori,itori1,itori2,iblock)
 +          v1cij=v1c(1,j,itori,itori1,itori2,iblock)
 +          v1sij=v1s(1,j,itori,itori1,itori2,iblock)
 +          v2cij=v1c(2,j,itori,itori1,itori2,iblock)
 +          v2sij=v1s(2,j,itori,itori1,itori2,iblock)
 +
            cosphi1=dcos(j*phii)
            sinphi1=dsin(j*phii)
            cosphi2=dcos(j*phii1)
@@@ -4558,12 -4512,12 +4561,12 @@@
            gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
            gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
          enddo
 -        do k=2,ntermd_2(itori,itori1,itori2)
 +        do k=2,ntermd_2(itori,itori1,itori2,iblock)
            do l=1,k-1
 -            v1cdij = v2c(k,l,itori,itori1,itori2)
 -            v2cdij = v2c(l,k,itori,itori1,itori2)
 -            v1sdij = v2s(k,l,itori,itori1,itori2)
 -            v2sdij = v2s(l,k,itori,itori1,itori2)
 +            v1cdij = v2c(k,l,itori,itori1,itori2,iblock)
 +            v2cdij = v2c(l,k,itori,itori1,itori2,iblock)
 +            v1sdij = v2s(k,l,itori,itori1,itori2,iblock)
 +            v2sdij = v2s(l,k,itori,itori1,itori2,iblock)
              cosphi1p2=dcos(l*phii+(k-l)*phii1)
              cosphi1m2=dcos(l*phii-(k-l)*phii1)
              sinphi1p2=dsin(l*phii+(k-l)*phii1)
@@@ -4614,8 -4568,8 +4617,8 @@@ c      write (iout,*) "EBACK_SC_COR",it
        esccor=0.0D0
        do i=itau_start,itau_end
          esccor_ii=0.0D0
 -        isccori=isccortyp(itype(i-2))
 -        isccori1=isccortyp(itype(i-1))
 +        isccori=isccortyp((itype(i-2)))
 +        isccori1=isccortyp((itype(i-1)))
          phii=phi(i)
  cccc  Added 9 May 2012
  cc Tauangle is torsional engle depending on the value of first digit 
@@@ -4632,14 -4586,14 +4635,14 @@@ c   2 = Ca...Ca...Ca...S
  c   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.21).or.
 -     &      (itype(i-1).eq.21)))
 +     &      (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.21)))
 +     &     .or.(itype(i-2).eq.ntyp1)))
       &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
 -     &      (itype(i-1).eq.21)))) cycle
 -        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle
 -        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21))
 +     &      (itype(i-1).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))
       & cycle
          do j=1,nterm_sccor(isccori,isccori1)
            v1ij=v1sccor(j,intertyp,isccori,isccori1)
@@@ -4657,7 -4611,7 +4660,7 @@@ c     &gloc_sc(intertyp,i-3,icg
       &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
       &  (v1sccor(j,intertyp,itori,itori1),j=1,6)
       & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
 -        gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
 +c        gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
         enddo !intertyp
        enddo
  c        do i=1,nres
@@@ -4770,9 -4724,9 +4773,9 @@@ c--------------------------------------
        integer dimen1,dimen2,atom,indx
        double precision buffer(dimen1,dimen2)
        double precision zapas 
 -      common /contacts_hb/ zapas(3,20,maxres,7),
 -     &   facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
 -     &         num_cont_hb(maxres),jcont_hb(20,maxres)
 +      common /contacts_hb/ zapas(3,ntyp,maxres,7),
 +     &   facont_hb(ntyp,maxres),ees0p(ntyp,maxres),ees0m(ntyp,maxres),
 +     &         num_cont_hb(maxres),jcont_hb(ntyp,maxres)
        num_kont=num_cont_hb(atom)
        do i=1,num_kont
          do k=1,7
@@@ -4795,10 -4749,9 +4798,10 @@@ c--------------------------------------
        integer dimen1,dimen2,atom,indx
        double precision buffer(dimen1,dimen2)
        double precision zapas 
 -      common /contacts_hb/ zapas(3,20,maxres,7),
 -     &         facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
 -     &         num_cont_hb(maxres),jcont_hb(20,maxres)
 +      common /contacts_hb/ zapas(3,ntyp,maxres,7),
 +     &         facont_hb(ntyp,maxres),ees0p(ntyp,maxres),
 +     &         ees0m(ntyp,maxres),
 +     &         num_cont_hb(maxres),jcont_hb(ntyp,maxres)
        num_kont=buffer(1,indx+26)
        num_kont_old=num_cont_hb(atom)
        num_cont_hb(atom)=num_kont+num_kont_old