Merge branch 'prerelease-3.2.1' into devel
authorDawid Jagiela <lightnir@chem.univ.gda.pl>
Wed, 5 Nov 2014 22:44:32 +0000 (23:44 +0100)
committerDawid Jagiela <lightnir@chem.univ.gda.pl>
Wed, 5 Nov 2014 22:44:32 +0000 (23:44 +0100)
Conflicts:
.gitignore
CMakeLists.txt
source/cluster/wham/src/COMMON.SCCOR
source/maxlik/src_CSA/CMakeLists.txt
source/unres/src_CSA_DiL/CMakeLists.txt
source/unres/src_MD-M/CMakeLists.txt
source/unres/src_MD-M/Makefile
source/unres/src_MD-M/energy_p_new_barrier.F
source/unres/src_MD-M/ssMD.F
source/unres/src_MD/CMakeLists.txt
source/unres/src_MIN/CMakeLists.txt
source/wham/src/DIMENSIONS.FREE
source/xdrfpdb/src/CMakeLists.txt

13 files changed:
1  2 
CMakeLists.txt
source/maxlik/src_CSA/CMakeLists.txt
source/unres/src_MD-M/CMakeLists.txt
source/unres/src_MD-M/energy_p_new_barrier.F
source/unres/src_MD-M/readrtns_CSA.F
source/unres/src_MD/CMakeLists.txt
source/unres/src_MD/COMMON.SCCOR
source/unres/src_MD/energy_p_new_barrier.F
source/unres/src_MIN/CMakeLists.txt
source/wham/src-M/CMakeLists.txt
source/wham/src/CMakeLists.txt
source/wham/src/energy_p_new.F
source/xdrfpdb/src/CMakeLists.txt

diff --combined 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 
  #=======================================
@@@ -125,6 -143,7 +143,7 @@@ find_package(MPI QUIET
  
  if (MPI_Fortran_FOUND)
    message("MPI found")
+   FIX_DBL_INCLUDE(MPI_Fortran_INCLUDE_PATH)
  else()
    message("MPI not found - disabling MPI compile flags ")
    set ( UNRES_WITH_MPI "OFF")
@@@ -149,26 -168,20 +168,19 @@@ 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)
@@@ -190,7 -203,6 +202,7 @@@ 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)
@@@ -38,7 -38,7 +38,7 @@@ endif(${CMAKE_SYSTEM_NAME} MATCHES "Lin
  #=========================================
  # Set binary name 
  #=========================================
- set(MAXLIK_BIN "maxlik-opt-multprot")
+ set(MAXLIK_BIN "maxlik_CSA")
  
  
  #=========================================
@@@ -52,9 -52,10 +52,11 @@@ set(MAXLIK_SRCS ${MAXLIK_SRC0} 
  #=========================================
  add_executable(MAXLIK ${MAXLIK_SRCS} )
  set_target_properties(MAXLIK PROPERTIES OUTPUT_NAME ${MAXLIK_BIN})
+ set_property(TARGET MAXLIK PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  
  #=========================================
 -# Install Path
 +# Install maxlik binary (used by "make install")
  #=========================================
+ install(TARGETS MAXLIK DESTINATION ${CMAKE_INSTALL_PREFIX})
  
 +install(TARGETS MAXLIK RUNTIME DESTINATION unrespack/bin)
@@@ -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 
        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
        newconf.f 
        parmread.F 
        permut.F
+       prng_32.F
        q_measure1.F 
        q_measure3.F 
        q_measure.F
@@@ -204,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 -DTIMING -DTIMING_ENE" )
  
  #=========================================
  #  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")
  
  
@@@ -236,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")
  
  
@@@ -250,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")
@@@ -269,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_MD_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_MD_FF}.exe")
+   set(UNRES_BIN "unresMD-M_${Fortran_COMPILER_NAME}_single_${UNRES_MD_FF}.exe")
  endif(UNRES_WITH_MPI)  
  
  #=========================================
@@@ -317,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})
  
  #=========================================
@@@ -333,10 -319,9 +322,9 @@@ endif(UNRES_WITH_MPI
  target_link_libraries( UNRES_BIN-MD-M xdrf )
  
  #=========================================
- # INSTALL
+ # Install Path
  #=========================================
- install(TARGETS UNRES_BIN-MD-M RUNTIME DESTINATION unrespack/bin)
+ install(TARGETS UNRES_BIN-MD-M DESTINATION ${CMAKE_INSTALL_PREFIX})
  
  #=========================================
  # TESTS 
        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
@@@ -442,7 -440,7 +442,7 @@@ cMS$ATTRIBUTES C ::  proc_pro
  #ifdef MPI
        include 'mpif.h'
        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.CONTROL'
        include 'COMMON.TIME1'
        include 'COMMON.MAXGRAD'
 +      include 'COMMON.SCCOR'
  #ifdef TIMING
        time01=MPI_Wtime()
  #endif
@@@ -698,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"
          do i=1,4*nres
            glocbuf(i)=gloc(i,icg)
          enddo
 +c#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
 +c#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
          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
 +c#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
 +c#undef DEBUG
  #ifdef DEBUG
        write (iout,*) "gloc after reduce"
        do i=1,4*nres
  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)
  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
  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)
@@@ -1225,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
@@@ -1307,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)
          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)
@@@ -1428,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)
       &                        '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,
@@@ -1580,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)
          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)
  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)
  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
@@@ -1910,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)
          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)
@@@ -2385,11 -2355,7 +2385,11 @@@ c        if (i .gt. iatel_s+2) the
          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
  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)
          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)
          zmedi=c(3,i)+0.5d0*dzi
          num_conti=num_cont_hb(i)
          call eelecij(i,i+3,ees,evdw1,eel_loc)
 -        if (wturn4.gt.0.0d0 .and. itype(i+2).ne.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
  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
 +        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
          dxi=dc(1,i)
          dyi=dc(2,i)
          dzi=dc(3,i)
@@@ -2856,7 -2822,7 +2856,7 @@@ c        write (iout,*) 'i',i,' ielstar
          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
 +          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
@@@ -2958,9 -2924,7 +2958,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
  
@@@ -3292,7 -3256,6 +3292,7 @@@ cd          write (iout,*) 'i',i,' j',j
  
            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
  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))
          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
  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))
          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
@@@ -3979,9 -3942,8 +3979,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
@@@ -4081,48 -4043,48 +4081,54 @@@ cmc        if (ii.gt.nres .and. itype(i
  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
++<<<<<<< HEAD
 +         if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
 +     & iabs(itype(jjj)).eq.1) then
++=======
+          if (ii.gt.nres 
+      &       .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then 
++>>>>>>> prerelease-3.2.1
            call ssbond_ene(iii,jjj,eij)
            ehpb=ehpb+2*eij
+          endif
  cd          write (iout,*) "eij",eij
 +         endif
          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
        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)
        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
        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)
  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)
@@@ -4336,24 -4298,11 +4342,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
            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
@@@ -4386,27 -4335,15 +4392,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)
        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)
              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
  #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)
              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)
          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
          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)
          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)
            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
@@@ -4770,9 -4699,9 +4776,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
@@@ -4929,11 -4858,11 +4935,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
@@@ -5015,7 -4944,7 +5021,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
        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)))
          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
@@@ -5096,7 -5025,7 +5102,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)
  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
  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
@@@ -5178,9 -5107,7 +5184,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
@@@ -5224,7 -5151,6 +5230,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
       &        +(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
       &        +(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
       &  +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))
@@@ -5310,16 -5235,13 +5316,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)
@@@ -5504,8 -5426,8 +5510,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)
@@@ -5601,22 -5523,17 +5607,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
@@@ -5631,7 -5548,7 +5637,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)
            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
@@@ -5701,10 -5617,9 +5707,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))
          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)
            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)
@@@ -5780,53 -5692,29 +5786,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----------------------------------------------------------------------------
  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, 
@@@ -8216,10 -8104,10 +8222,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
@@@ -8333,7 -8221,7 +8339,7 @@@ c--------------------------------------
       & auxvec1(2),auxmat1(2,2)
        logical swap
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- C                                                                              C
+ C                                                                              C                       
  C      Parallel       Antiparallel                                             C
  C                                                                              C
  C          o             o                                                     C
@@@ -8341,10 -8229,10 +8347,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 
@@@ -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--------------------------------------------------------------------------
        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) 
  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,
       &  '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)
@@@ -741,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)
@@@ -773,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
@@@ -783,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
@@@ -848,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
@@@ -966,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)
           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)
@@@ -1088,8 -1088,7 +1089,7 @@@ C Generate distance constraints, if th
          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
@@@ -1165,7 -1164,7 +1165,7 @@@ 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.
          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)
@@@ -1337,7 -1336,7 +1337,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
@@@ -1961,22 -1960,12 +1961,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
@@@ -2138,8 -2127,6 +2138,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) 
@@@ -2281,11 -2268,11 +2281,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
  #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.
@@@ -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 
        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
        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,11 -145,6 +145,11 @@@ 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")
  
  
@@@ -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} )  
  #========================================
  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})
  
  
@@@ -308,10 -288,9 +299,9 @@@ endif(UNRES_WITH_MPI
  target_link_libraries( UNRES_BIN-MD xdrf )
  
  #=========================================
- # INSTALL 
+ # Install Path
  #=========================================
-  
- install(TARGETS UNRES_BIN-MD RUNTIME DESTINATION unrespack/bin)
+ install(TARGETS UNRES_BIN-MD DESTINATION ${CMAKE_INSTALL_PREFIX}) 
  
  #=========================================
  # TESTS 
  #=========================================
  #  test_single_ala.sh
  #=========================================
 +#
 +# Set parmaeters depending on force field
 +if(UNRES_MD_FF STREQUAL "GAB")
 +   set(UNRES_BONDPAR "bond.parm")     
 +elseif(UNRES_MD_FF STREQUAL "E0LL2Y")
 +   set(UNRES_BONDPAR "bond_AM1.parm")
 +endif(UNRES_MD_FF STREQUAL "GAB")
  
  FILE(WRITE ${CMAKE_CURRENT_BINARY_DIR}/scripts/test_single_ala.sh
  "#!/bin/sh
  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
 +export BONDPAR=$DD/${UNRES_BONDPAR}
  export THETPAR=$DD/thetaml.5parm
  export ROTPAR=$DD/scgauss.parm
  export TORPAR=$DD/torsion_631Gdp.parm
@@@ -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)
          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)
            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
              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)
          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
          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)
          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)
            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
          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)
  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
  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
@@@ -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,47 -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 )
  
  #=========================================
- # INSTALL 
+ # Install Path
  #=========================================
+ install(TARGETS UNRES_MIN_BIN DESTINATION ${CMAKE_INSTALL_PREFIX})
  
- install(TARGETS UNRES_BIN-MIN RUNTIME DESTINATION unrespack/bin)
- #=========================================
- # TESTS 
- #=========================================
- #-- 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)
-  
@@@ -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")
  
  
@@@ -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,8 -220,7 +224,7 @@@ 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})
  
  #=========================================
@@@ -234,6 -231,13 +235,13 @@@ target_link_libraries( UNRES_WHAM_M_BI
  # link libxdrf.a 
  target_link_libraries( UNRES_WHAM_M_BIN xdrf )
  
+ #=========================================
+ # Install Path
+ #=========================================
+ install(TARGETS UNRES_WHAM_M_BIN DESTINATION ${CMAKE_INSTALL_PREFIX})
  #=========================================
  # TESTS 
  #=========================================
@@@ -94,9 -94,11 +94,11 @@@ set(UNRES_WHAM_PP_SR
  # Set comipiler flags for different sourcefiles  
  #================================================
  if (Fortran_COMPILER_NAME STREQUAL "ifort")
 -  set(FFLAGS0 "-mcmodel=medium -g -CB -I. -I${CMAKE_CURRENT_SOURCE_DIR}/include_unres" ) 
 +  set(FFLAGS0 "-mcmodel=large -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")
  
  
@@@ -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,8 -221,7 +224,7 @@@ 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})
  
  #=========================================
@@@ -235,6 -233,12 +236,12 @@@ target_link_libraries( UNRES_WHAM_BIN $
  target_link_libraries( UNRES_WHAM_BIN xdrf )
  
  #=========================================
+ # Install Path
+ #=========================================
+ install(TARGETS UNRES_WHAM_BIN DESTINATION ${CMAKE_INSTALL_PREFIX})
+ #=========================================
  # TESTS 
  #=========================================
  
       &   +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)
  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)
          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)
          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
        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)
        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
  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
          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)
        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)
  #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)
  #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)
  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)
          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
          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)
          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)
            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
          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)
  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
  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)
            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)
            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
@@@ -71,31 -71,36 +71,36 @@@ set_property(SOURCE ${UNRES_XDRF_PP_SRC
  
   
  #=========================================
 -# Build the binaries
 +# Build the binaries (used by "make")
  #=========================================
  add_executable(xdrf2pdb   ${UNRES_XDRF_XDRF2PDB_SRC} )
+ set_property(TARGET xdrf2pdb PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  target_link_libraries( xdrf2pdb xdrf )
  
  add_executable(xdrf2pdb-m ${UNRES_XDRF_XDRF2PDB-M_SRC} )
+ set_property(TARGET xdrf2pdb-m PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  target_link_libraries( xdrf2pdb-m xdrf )
  
  add_executable(xdrf2x     ${UNRES_XDRF_XDRF2X_SRC} )
+ set_property(TARGET xdrf2x PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  target_link_libraries( xdrf2x xdrf )
  
  add_executable(xdrf2ang   ${UNRES_XDRF_XDRF2XANG_SRC} )
+ set_property(TARGET xdrf2ang PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
  target_link_libraries( xdrf2ang xdrf )
  
  
  #=========================================
- # INSTALL (used by "make install")
+ # Install Path
  #=========================================
- install(TARGETS xdrf2pdb      RUNTIME DESTINATION unrespack/bin)
- install(TARGETS xdrf2pdb-m    RUNTIME DESTINATION unrespack/bin)
- install(TARGETS xdrf2x                RUNTIME DESTINATION unrespack/bin)
- install(TARGETS xdrf2ang      RUNTIME DESTINATION unrespack/bin)
 -install(TARGETS xdrf2pdb DESTINATION ${CMAKE_INSTALL_PREFIX})
 -install(TARGETS xdrf2pdb-m DESTINATION ${CMAKE_INSTALL_PREFIX})
 -install(TARGETS xdrf2x DESTINATION ${CMAKE_INSTALL_PREFIX})
 -install(TARGETS xdrf2ang DESTINATION ${CMAKE_INSTALL_PREFIX})
++install(TARGETS xdrf2pdb      DESTINATION ${CMAKE_INSTALL_PREFIX})
++install(TARGETS xdrf2pdb-m    DESTINATION ${CMAKE_INSTALL_PREFIX})
++install(TARGETS xdrf2x                DESTINATION ${CMAKE_INSTALL_PREFIX})
++install(TARGETS xdrf2ang      DESTINATION ${CMAKE_INSTALL_PREFIX})
  
  #=========================================
- # TESTS (used by "ctest") 
+ # TESTS
  #=========================================