Adam's changes
[unres.git] / source / unres / src-HCD-5D / readrtns_CSA.F
index 66f7c17..d76b29e 100644 (file)
@@ -1,5 +1,5 @@
       subroutine readrtns
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
@@ -9,6 +9,7 @@
       include 'COMMON.SBRIDGE'
       include 'COMMON.IOUNITS'
       include 'COMMON.SPLITELE'
+      integer i,j
       logical file_exist
 C Read job setup parameters
       call read_control
@@ -84,18 +85,19 @@ C-------------------------------------------------------------------------------
 C
 C Read contorl data
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MP
       include 'mpif.h'
       logical OKRandom, prng_restart
-      real*8  r1
+      double precision  r1
 #endif
       include 'COMMON.IOUNITS'
       include 'COMMON.TIME1'
       include 'COMMON.THREAD'
       include 'COMMON.SBRIDGE'
       include 'COMMON.CONTROL'
+      include 'COMMON.SAXS'
       include 'COMMON.MCM'
       include 'COMMON.MAP'
       include 'COMMON.HEADER'
@@ -109,10 +111,13 @@ C
       include 'COMMON.SPLITELE'
       include 'COMMON.SHIELD'
       include 'COMMON.GEO'
+      integer i
+      integer KDIAG,ICORFL,IXDR
       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
       character*80 ucase
       character*320 controlcard
+      double precision seed
 
       nglob_csa=0
       eglob_csa=1d99
@@ -125,6 +130,9 @@ c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
       call random_init(seed)
 C Set up the time limit (caution! The time must be input in minutes!)
       read_cart=index(controlcard,'READ_CART').gt.0
+      out_cart=index(controlcard,'OUT_CART').gt.0
+      out_int=index(controlcard,'OUT_INT').gt.0
+      gmatout=index(controlcard,'GMATOUT').gt.0
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 C this variable with_theta_constr is the variable which allow to read and execute the
 C constrains on theta angles WITH_THETA_CONSTR is the keyword
@@ -144,33 +152,33 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword
       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
-      call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
-      call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
-      call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
-      call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
-      call reada(controlcard,'DRMS',drms,0.1D0)
-      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
-       write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc 
-       write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1 
-       write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max 
-       write (iout,'(a,f10.1)')'DRMS    = ',drms 
-       write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm 
-       write (iout,'(a,f10.1)') 'Time limit (min):',timlim
-      endif
+c      call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
+c      call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
+c      call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
+c      call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
+c      call reada(controlcard,'DRMS',drms,0.1D0)
+c      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+c       write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc 
+c       write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1 
+c       write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max 
+c       write (iout,'(a,f10.1)')'DRMS    = ',drms 
+cc       write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm 
+c       write (iout,'(a,f10.1)') 'Time limit (min):',timlim
+c      endif
       call readi(controlcard,'NZ_START',nz_start,0)
       call readi(controlcard,'NZ_END',nz_end,0)
 c      call readi(controlcard,'IZ_SC',iz_sc,0)
       timlim=60.0D0*timlim
       safety = 60.0d0*safety
-      timem=timlim
       modecalc=0
+      call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
       call reada(controlcard,"T_BATH",t_bath,300.0d0)
       minim=(index(controlcard,'MINIMIZE').gt.0)
       dccart=(index(controlcard,'CART').gt.0)
       overlapsc=(index(controlcard,'OVERLAP').gt.0)
       overlapsc=.not.overlapsc
-      searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
-      searchsc=.not.searchsc
+      searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
+     &  index(controlcard,'NOSEARCHSC').eq.0)
       sideadd=(index(controlcard,'SIDEADD').gt.0)
       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
       mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
@@ -179,7 +187,6 @@ c      call readi(controlcard,'IZ_SC',iz_sc,0)
       pdbref=(index(controlcard,'PDBREF').gt.0)
       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
       indpdb=index(controlcard,'PDBSTART')
-      extconf=(index(controlcard,'EXTCONF').gt.0)
       AFMlog=(index(controlcard,'AFM'))
       selfguide=(index(controlcard,'SELFGUIDE'))
 c      print *,'AFMlog',AFMlog,selfguide,"KUPA"
@@ -288,6 +295,12 @@ cfmc        modecalc=10
       indphi=index(controlcard,'PHI')
       indback=index(controlcard,'BACK')
       iranconf=index(controlcard,'RAND_CONF')
+      start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
+      extconf=(index(controlcard,'EXTCONF').gt.0)
+      if (start_from_model) then
+        iranconf=0
+        extconf=.false.
+      endif
       i2ndstr=index(controlcard,'USE_SEC_PRED')
       gradout=index(controlcard,'GRADOUT').gt.0
       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
@@ -297,9 +310,16 @@ C Reading the dimensions of box in x,y,z coordinates
       call reada(controlcard,'BOXX',boxxsize,100.0d0)
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+      write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
 c Cutoff range for interactions
-      call reada(controlcard,"R_CUT",r_cut,15.0d0)
+      call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
+      call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+      write (iout,*) "Cutoff on interactions",r_cut_int
+      write (iout,*) 
+     & "Cutoff in switching short and long range interactions in RESPA",
+     & r_cut_respa
+      write (iout,*) "lambda in switch function",rlamb
       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
       if (lipthick.gt.0.0d0) then
@@ -328,25 +348,6 @@ C      endif
        buftubebot=bordtubebot+tubebufthick
        buftubetop=bordtubetop-tubebufthick
       endif
-c      if (shield_mode.gt.0) then
-c      pi=3.141592d0
-C VSolvSphere the volume of solving sphere
-C      print *,pi,"pi"
-C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
-C there will be no distinction between proline peptide group and normal peptide
-C group in case of shielding parameters
-c      write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
-c      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
-c      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
-c      write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
-c     &  VSolvSphere_div
-C long axis of side chain 
-c      do i=1,ntyp
-c      long_r_sidechain(i)=vbldsc0(1,i)
-c      short_r_sidechain(i)=sigma0(i)
-c      enddo
-c      buff_shield=1.0d0
-c      endif
       if (me.eq.king .or. .not.out1file ) 
      &  write (iout,*) "DISTCHAINMAX",distchainmax
       
@@ -360,7 +361,7 @@ c--------------------------------------------------------------------------
 C
 C Read REMD settings
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.IOUNITS'
       include 'COMMON.TIME1'
@@ -368,8 +369,12 @@ C
 #ifndef LANG0
       include 'COMMON.LANGEVIN'
 #else
+#ifdef FIVEDIAG
+      include 'COMMON.LANGEVIN.lang0.5diag'
+#else
       include 'COMMON.LANGEVIN.lang0'
 #endif
+#endif
       include 'COMMON.INTERACT'
       include 'COMMON.NAMES'
       include 'COMMON.GEO'
@@ -379,7 +384,7 @@ C
       character*80 ucase
       character*320 controlcard
       character*3200 controlcard1
-      integer iremd_m_total
+      integer iremd_m_total,i
 
       if(me.eq.king.or..not.out1file)
      & write (iout,*) "REMD setup"
@@ -445,16 +450,21 @@ c--------------------------------------------------------------------------
 C
 C Read MD settings
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.IOUNITS'
       include 'COMMON.TIME1'
       include 'COMMON.MD'
+      include 'COMMON.QRESTR'
 #ifndef LANG0
       include 'COMMON.LANGEVIN'
 #else
+#ifdef FIVEDIAG
+      include 'COMMON.LANGEVIN.lang0.5diag'
+#else
       include 'COMMON.LANGEVIN.lang0'
 #endif
+#endif
       include 'COMMON.INTERACT'
       include 'COMMON.NAMES'
       include 'COMMON.GEO'
@@ -464,6 +474,8 @@ C
       include 'COMMON.FFIELD'
       character*80 ucase
       character*320 controlcard
+      integer i
+      double precision eta
 
       call card_concat(controlcard)
       call readi(controlcard,"NSTEP",n_timestep,1000000)
@@ -540,12 +552,16 @@ c  if performing umbrella sampling, fragments constrained are read from the frag
      &  "Initial time step of numerical integration:",d_time,
      &  " natural units"
        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
+       write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
+     &   " A"
+       write(iout,'(a60,i5)')"Frequency of updating interaction list",
+     &   imatupdate
        if (RESPA) then
         write (iout,'(2a,i4,a)') 
      &    "A-MTS algorithm used; initial time step for fast-varying",
      &    " short-range forces split into",ntime_split," steps."
         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
-     &   r_cut," lambda",rlamb
+     &   r_cut_respa," lambda",rlamb
        endif
        write (iout,'(2a,f10.5)') 
      &  "Maximum acceleration threshold to reduce the time step",
@@ -681,11 +697,11 @@ c------------------------------------------------------------------------------
 C
 C Read molecular data.
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
-      integer error_msg
+      integer error_msg,ierror,ierr,ierrcode
 #endif
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
@@ -698,6 +714,7 @@ C
       include 'COMMON.SBRIDGE'
       include 'COMMON.HEADER'
       include 'COMMON.CONTROL'
+      include 'COMMON.SAXS'
       include 'COMMON.DBASE'
       include 'COMMON.THREAD'
       include 'COMMON.CONTACTS'
@@ -713,14 +730,19 @@ C
       character*256 pdbfile
       character*400 weightcard
       character*80 weightcard_t,ucase
-      dimension itype_pdb(maxres)
+      integer itype_pdb(maxres)
       common /pizda/ itype_pdb
       logical seq_comp,fail
       double precision energia(0:n_ene)
       double precision secprob(3,maxdih_constr)
+      double precision co
+      double precision phihel,phibet,sigmahel,sigmabet
+      integer iti,nsi,maxsi
       integer ilen
       external ilen
-      integer tperm
+      integer iperm,tperm
+      integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
+      double precision sumv
 C
 C Read PDB structure if applicable
 C
@@ -789,25 +811,6 @@ C Convert sequence to numeric code
         do i=1,nres
           itype(i)=rescode(i,sequence(i),iscode)
         enddo
-C Assign initial virtual bond lengths
-c        do i=2,nres
-c          vbld(i)=vbl
-c          vbld_inv(i)=vblinv
-c        enddo
-c        if (itype(1).eq.ntyp1) then
-c          vbld(2)=vbld(2)/2
-c          vbld_inv(2)=vbld_inv(2)*2
-c        endif
-c        if (itype(nres).eq.ntyp1) then
-c          vbld(nres)=vbld(nres)/2
-c          vbld_inv(nres)=vbld_inv(nres)*2
-c        endif
-c        do i=2,nres-1
-c          vbld(i+nres)=dsc(iabs(itype(i)))
-c          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)
-c        enddo
       endif 
 c      print *,nres
 c      print '(20i4)',(itype(i),i=1,nres)
@@ -833,17 +836,25 @@ c      print '(20i4)',(itype(i),i=1,nres)
        do i=1,nres-1
          write (iout,*) i,itype(i),itel(i)
        enddo
-       print *,'Call Read_Bridge.'
+c       print *,'Call Read_Bridge.'
       endif
       nnt=1
       nct=nres
 cd      print *,'NNT=',NNT,' NCT=',NCT
       call seq2chains(nres,itype,nchain,chain_length,chain_border,
      &  ireschain)
+      chain_border1(1,1)=1
+      chain_border1(2,1)=chain_border(2,1)+1
+      do i=2,nchain-1
+        chain_border1(1,i)=chain_border(1,i)-1
+        chain_border1(2,i)=chain_border(2,i)+1
+      enddo
+      if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
+      chain_border1(2,nchain)=nres
       write(iout,*) "nres",nres," nchain",nchain
       do i=1,nchain
         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
-     &    chain_border(2,i)
+     &    chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
       enddo
       call chain_symmetry(nchain,nres,itype,chain_border,
      &    chain_length,npermchain,tabpermchain)
@@ -855,15 +866,18 @@ c      enddo
       do i=1,nres
         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
       enddo
+      call flush(iout)
       if (itype(1).eq.ntyp1) nnt=2
       if (itype(nres).eq.ntyp1) nct=nct-1
+      write (iout,*) "nnt",nnt," nct",nct
+      call flush(iout)
 #ifdef DFA
       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
        call init_dfa_vars
-       print*, 'init_dfa_vars finished!'
+c       print*, 'init_dfa_vars finished!'
        call read_dfa_info
-       print*, 'read_dfa_info finished!'
+c       print*, 'read_dfa_info finished!'
       endif
 #endif
       if (pdbref) then
@@ -1088,10 +1102,10 @@ czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
         endif
       endif
 c        print *, "A TU"
-      write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
+c      write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
       call flush(iout)
       if (constr_dist.gt.0) call read_dist_constr
-      write (iout,*) "After read_dist_constr nhpb",nhpb
+c      write (iout,*) "After read_dist_constr nhpb",nhpb
       if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
       call hpb_partition
       call NMRpeak_partition
@@ -1154,6 +1168,49 @@ c        print *, "A TU"
        enddo
       else
         homol_nset=0
+        if (start_from_model) then
+          nmodel_start=0 
+          do 
+            read(inp,'(a)',end=332,err=332) pdbfile
+            if (me.eq.king .or. .not. out1file)
+     &        write (iout,'(a,5x,a)') 'Opening PDB file',
+     &        pdbfile(:ilen(pdbfile))
+            open(ipdbin,file=pdbfile,status='old',err=336)
+            goto 335
+ 336        write (iout,'(a,5x,a)') 'Error opening PDB file',
+     &      pdbfile(:ilen(pdbfile))
+            call flush(iout)
+            stop
+ 335        continue
+            unres_pdb=.false.
+            nres_temp=nres
+            call readpdb
+            close(ipdbin)
+            if (nres.ge.nres_temp) then
+              nmodel_start=nmodel_start+1
+              pdbfiles_chomo(nmodel_start)=pdbfile
+              do i=1,2*nres
+                do j=1,3
+                  chomo(j,i,nmodel_start)=c(j,i)
+                enddo
+              enddo
+            else
+              if (me.eq.king .or. .not. out1file)
+     &          write (iout,'(a,2i5,1x,a)') 
+     &           "Different number of residues",nres_temp,nres,
+     &           " model skipped."
+            endif
+            nres=nres_temp
+          enddo
+  332     continue
+          if (nmodel_start.eq.0) then
+            if (me.eq.king .or. .not. out1file)
+     &        write (iout,'(a)')
+     &        "No valid starting model found START_FROM_MODELS is OFF"
+              start_from_model=.false.
+          endif
+          write (iout,*) "nmodel_start",nmodel_start
+        endif
       endif
 
 
@@ -1163,14 +1220,15 @@ C      endif
      &    modecalc.ne.10) then
 C If input structure hasn't been supplied from the PDB file read or generate
 C initial geometry.
-        if (iranconf.eq.0 .and. .not. extconf) then
+        if (iranconf.eq.0 .and. .not. extconf .and. .not.
+     &     start_from_model) then
           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
      &     write (iout,'(a)') 'Initial geometry will be read in.'
           if (read_cart) then
             read(inp,'(8f10.5)',end=36,err=36)
      &       ((c(l,k),l=1,3),k=1,nres),
      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
-            write (iout,*) "Exit READ_CART"
+c            write (iout,*) "Exit READ_CART"
 c            write (iout,'(8f10.5)') 
 c     &       ((c(l,k),l=1,3),k=1,nres),
 c     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
@@ -1372,10 +1430,11 @@ c--------------------------------------------------------------------------
 c-----------------------------------------------------------------------------
       subroutine read_bridge
 C Read information about disulfide bridges.
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
+      integer ierror
 #endif
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
@@ -1392,9 +1451,16 @@ C Read information about disulfide bridges.
       include 'COMMON.THREAD'
       include 'COMMON.TIME1'
       include 'COMMON.SETUP'
+      integer i,j
 C Read bridging residues.
       read (inp,*) ns,(iss(i),i=1,ns)
-      print *,'ns=',ns
+c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
+c     numbers
+      icys=0
+      do i=1,ns
+        icys(iss(i))=i
+      enddo
+c      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.
@@ -1447,8 +1513,8 @@ C bridging residues.
           enddo
           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
    20     continue
-          dhpb(i)=dbr
-          forcon(i)=fbr
+c          dhpb(i)=dbr
+c          forcon(i)=fbr
         enddo
         do i=1,nss
           ihpb(i)=ihpb(i)+nres
@@ -1460,7 +1526,7 @@ C bridging residues.
       end
 c----------------------------------------------------------------------------
       subroutine read_x(kanal,*)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -1469,6 +1535,7 @@ c----------------------------------------------------------------------------
       include 'COMMON.CONTROL'
       include 'COMMON.LOCAL'
       include 'COMMON.INTERACT'
+      integer i,j,k,l,kanal
 c Read coordinates from input
 c
       read(kanal,'(8f10.5)',end=10,err=10)
@@ -1499,7 +1566,7 @@ c
       end
 c----------------------------------------------------------------------------
       subroutine read_threadbase
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
@@ -1515,6 +1582,8 @@ c----------------------------------------------------------------------------
       include 'COMMON.DBASE'
       include 'COMMON.THREAD'
       include 'COMMON.TIME1'
+      integer i,j,k
+      double precision dist
 C Read pattern database for threading.
       read (icbase,*) nseq
       do i=1,nseq
@@ -1544,7 +1613,7 @@ c    &   nres_base(1,i))
       end
 c------------------------------------------------------------------------------
       subroutine setup_var
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
@@ -1560,17 +1629,17 @@ c------------------------------------------------------------------------------
       include 'COMMON.DBASE'
       include 'COMMON.THREAD'
       include 'COMMON.TIME1'
+      integer i
 C Set up variable list.
       ntheta=nres-2
       nphi=nres-3
       nvar=ntheta+nphi
       nside=0
-      write (iout,*) "SETUP_VAR ialph"
       do i=2,nres-1
         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
-         nside=nside+1
+          nside=nside+1
           ialph(i,1)=nvar+nside
-         ialph(nside,2)=i
+          ialph(nside,2)=i
         endif
       enddo
       if (indphi.gt.0) then
@@ -1580,13 +1649,12 @@ C Set up variable list.
       else
         nvar=nvar+2*nside
       endif
-      write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
       return
       end
 c----------------------------------------------------------------------------
       subroutine gen_dist_constr
 C Generate CA distance constraints.
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
@@ -1601,9 +1669,12 @@ C Generate CA distance constraints.
       include 'COMMON.CONTROL'
       include 'COMMON.DBASE'
       include 'COMMON.THREAD'
+      include 'COMMON.SPLITELE'
       include 'COMMON.TIME1'
-      dimension itype_pdb(maxres)
+      integer i,j,itype_pdb(maxres)
       common /pizda/ itype_pdb
+      double precision dd
+      double precision dist
       character*2 iden
 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
@@ -1613,11 +1684,14 @@ cd     & ' nsup',nsup
 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
 cd     &    ' seq_pdb', restyp(itype_pdb(i))
         do j=i+2,nstart_sup+nsup-1
+c 5/24/2020 Adam: Cutoff included to reduce array size
+          dd = dist(i,j)
+          if (dd.gt.r_cut_int) cycle
           nhpb=nhpb+1
           ihpb(nhpb)=i+nstart_seq-nstart_sup
           jhpb(nhpb)=j+nstart_seq-nstart_sup
           forcon(nhpb)=weidis
-          dhpb(nhpb)=dist(i,j)
+          dhpb(nhpb)=dd
         enddo
       enddo 
 cd      write (iout,'(a)') 'Distance constraints:' 
@@ -1638,10 +1712,11 @@ cd      enddo
       end
 c----------------------------------------------------------------------------
       subroutine map_read
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.MAP'
       include 'COMMON.IOUNITS'
+      integer imap
       character*3 angid(4) /'THE','PHI','ALP','OME'/
       character*80 mapcard,ucase
       do imap=1,nmap
@@ -1684,7 +1759,7 @@ c----------------------------------------------------------------------------
       end 
 c----------------------------------------------------------------------------
       subroutine csaread
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
@@ -1746,114 +1821,15 @@ c!bankt
       return
       end
 c----------------------------------------------------------------------------
-cfmc      subroutine mcmfread
-cfmc      implicit real*8 (a-h,o-z)
-cfmc      include 'DIMENSIONS'
-cfmc      include 'COMMON.MCMF'
-cfmc      include 'COMMON.IOUNITS'
-cfmc      include 'COMMON.GEO'
-cfmc      character*80 ucase
-cfmc      character*620 mcmcard
-cfmc      call card_concat(mcmcard)
-cfmc
-cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
-cfmc      write(iout,*)'MAXRANT=',maxrant
-cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
-cfmc      write(iout,*)'MAXFAM=',maxfam
-cfmc      call readi(mcmcard,'NNET1',nnet1,5)
-cfmc      write(iout,*)'NNET1=',nnet1
-cfmc      call readi(mcmcard,'NNET2',nnet2,4)
-cfmc      write(iout,*)'NNET2=',nnet2
-cfmc      call readi(mcmcard,'NNET3',nnet3,4)
-cfmc      write(iout,*)'NNET3=',nnet3
-cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
-cfmc      write(iout,*)'ILASTT=',ilastt
-cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
-cfmc      write(iout,*)'MAXSTR=',maxstr
-cfmc      maxstr_f=maxstr/maxfam
-cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
-cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
-cfmc      write(iout,*)'NMCMF=',nmcmf
-cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
-cfmc      write(iout,*)'IFOCUS=',ifocus
-cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
-cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
-cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
-cfmc      write(iout,*)'INTPRT=',intprt
-cfmc      call readi(mcmcard,'IPRT',iprt,100)
-cfmc      write(iout,*)'IPRT=',iprt
-cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
-cfmc      write(iout,*)'IMAXTR=',imaxtr
-cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
-cfmc      write(iout,*)'MAXEVEN=',maxeven
-cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
-cfmc      write(iout,*)'MAXEVEN1=',maxeven1
-cfmc      call readi(mcmcard,'INIMIN',inimin,200)
-cfmc      write(iout,*)'INIMIN=',inimin
-cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
-cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
-cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
-cfmc      write(iout,*)'NTHREAD=',nthread
-cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
-cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
-cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
-cfmc      write(iout,*)'MAXPERT=',maxpert
-cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
-cfmc      write(iout,*)'IRMSD=',irmsd
-cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
-cfmc      write(iout,*)'DENEMIN=',denemin
-cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
-cfmc      write(iout,*)'RCUT1S=',rcut1s
-cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
-cfmc      write(iout,*)'RCUT1E=',rcut1e
-cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
-cfmc      write(iout,*)'RCUT2S=',rcut2s
-cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
-cfmc      write(iout,*)'RCUT2E=',rcut2e
-cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
-cfmc      write(iout,*)'DPERT1=',d_pert1
-cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
-cfmc      write(iout,*)'DPERT1A=',d_pert1a
-cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
-cfmc      write(iout,*)'DPERT2=',d_pert2
-cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
-cfmc      write(iout,*)'DPERT2A=',d_pert2a
-cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
-cfmc      write(iout,*)'DPERT2B=',d_pert2b
-cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
-cfmc      write(iout,*)'DPERT2C=',d_pert2c
-cfmc      d_pert1=deg2rad*d_pert1
-cfmc      d_pert1a=deg2rad*d_pert1a
-cfmc      d_pert2=deg2rad*d_pert2
-cfmc      d_pert2a=deg2rad*d_pert2a
-cfmc      d_pert2b=deg2rad*d_pert2b
-cfmc      d_pert2c=deg2rad*d_pert2c
-cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
-cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
-cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
-cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
-cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
-cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
-cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
-cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
-cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
-cfmc      write(iout,*)'RCUTINI=',rcutini
-cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
-cfmc      write(iout,*)'GRAT=',grat
-cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
-cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
-cfmc
-cfmc      return
-cfmc      end 
-c----------------------------------------------------------------------------
       subroutine mcmread
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.MCM'
       include 'COMMON.MCE'
       include 'COMMON.IOUNITS'
       character*80 ucase
       character*320 mcmcard
+      integer i
       call card_concat(mcmcard)
       call readi(mcmcard,'MAXACC',maxacc,100)
       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
@@ -1913,7 +1889,7 @@ C Probabilities of different move types
       end 
 c----------------------------------------------------------------------------
       subroutine read_minim
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.MINIM'
       include 'COMMON.IOUNITS'
@@ -1939,13 +1915,14 @@ c----------------------------------------------------------------------------
       end
 c----------------------------------------------------------------------------
       subroutine read_angles(kanal,*)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
+      integer i,kanal
 c Read angles from input 
 c
        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
@@ -2039,10 +2016,11 @@ c----------------------------------------------------------------------------
       end
 c----------------------------------------------------------------------------
       subroutine openunits
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'    
 #ifdef MPI
       include 'mpif.h'
+      integer ierror
       character*16 form,nodename
       integer nodelen
 #endif
@@ -2050,7 +2028,7 @@ c----------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.MD'
       include 'COMMON.CONTROL'
-      integer lenpre,lenpot,ilen,lentmp
+      integer lenpre,lenpot,ilen,lentmp,npos
       external ilen
       character*3 out1file_text,ucase
       character*3 ll
@@ -2439,20 +2417,28 @@ c----------------------------------------------------------------------------
       card=card(:ilen(card)+1)//karta
       return
       end
-c----------------------------------------------------------------------------------
+c------------------------------------------------------------------------------
       subroutine readrst
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
+      include 'COMMON.CONTROL'
       include 'COMMON.MD'
+      include 'COMMON.QRESTR'
+#ifdef FIVEDIAG
+       include 'COMMON.LAGRANGE.5diag'
+#else
+       include 'COMMON.LAGRANGE'
+#endif
+      integer i,j
       open(irest2,file=rest2name,status='unknown')
       read(irest2,*) totT,EK,potE,totE,t_bath
       totTafm=totT
-      do i=1,2*nres
+      do i=0,2*nres-1
          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
       enddo
-      do i=1,2*nres
+      do i=0,2*nres-1
          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
       enddo
       if(usampl) then
@@ -2461,9 +2447,9 @@ c-------------------------------------------------------------------------------
       close(irest2)
       return
       end
-c---------------------------------------------------------------------------------
+c------------------------------------------------------------------------------
       subroutine read_fragments
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
@@ -2472,7 +2458,9 @@ c-------------------------------------------------------------------------------
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
       include 'COMMON.MD'
+      include 'COMMON.QRESTR'
       include 'COMMON.CONTROL'
+      integer i
       read(inp,*) nset,nfrag,npair,nfrag_back
       loc_qlike=(nfrag_back.lt.0)
       nfrag_back=iabs(nfrag_back)
@@ -2522,7 +2510,7 @@ c-------------------------------------------------------------------------------
       end
 C---------------------------------------------------------------------------
       subroutine read_afminp
-            implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
@@ -2533,35 +2521,38 @@ C---------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.SBRIDGE'
       character*320 afmcard
-      print *, "wchodze"
+      integer i
+c      print *, "wchodze"
       call card_concat(afmcard)
       call readi(afmcard,"BEG",afmbeg,0)
       call readi(afmcard,"END",afmend,0)
       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
       call reada(afmcard,"VEL",velAFMconst,0.0d0)
-      print *,'FORCE=' ,forceAFMconst
+c      print *,'FORCE=' ,forceAFMconst
 CCCC NOW PROPERTIES FOR AFM
        distafminit=0.0d0
        do i=1,3
         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
        enddo
         distafminit=dsqrt(distafminit)
-        print *,'initdist',distafminit
+c        print *,'initdist',distafminit
       return
       end
 c-------------------------------------------------------------------------------
       subroutine read_saxs_constr
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
 #endif
       include 'COMMON.SETUP'
       include 'COMMON.CONTROL'
+      include 'COMMON.SAXS'
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
       include 'COMMON.SBRIDGE'
-      double precision cm(3)
+      double precision cm(3),cnorm
+      integer i,j
 c      read(inp,*) nsaxs
       write (iout,*) "Calling read_saxs nsaxs",nsaxs
       call flush(iout)
@@ -2618,7 +2609,7 @@ c SAXS "spheres".
 
 c-------------------------------------------------------------------------------
       subroutine read_dist_constr
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
@@ -2629,8 +2620,10 @@ c-------------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.SBRIDGE'
       include 'COMMON.INTERACT'
-      integer ifrag_(2,100),ipair_(2,1000)
+      integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
+      integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
       double precision wfrag_(100),wpair_(1000)
+      double precision ddjk,dist,dist_cut,fordepthmax
       character*5000 controlcard
       logical normalize,next
       integer restr_type
@@ -3017,16 +3010,18 @@ C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
       end
 c-------------------------------------------------------------------------------
       subroutine read_constr_homology
-
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
 #endif
       include 'COMMON.SETUP'
       include 'COMMON.CONTROL'
+      include 'COMMON.HOMOLOGY'
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
       include 'COMMON.MD'
+      include 'COMMON.QRESTR'
       include 'COMMON.GEO'
       include 'COMMON.INTERACT'
       include 'COMMON.NAMES'
@@ -3043,10 +3038,12 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
       character*2 kic2
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
-      integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
+      integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
+     & ik,iistart,nres_temp
       integer ilen
       external ilen
-      logical liiflag
+      logical liiflag,lfirst
+      integer i01,i10
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
 c
@@ -3055,6 +3052,7 @@ c
       double precision, dimension (max_template,maxres) :: rescore
       double precision, dimension (max_template,maxres) :: rescore2
       double precision, dimension (max_template,maxres) :: rescore3
+      double precision distal
       character*24 pdbfile,tpl_k_rescore
 c -----------------------------------------------------------------
 c Reading multiple PDB ref structures and calculation of retraints
@@ -3078,15 +3076,17 @@ c Alternative: reading from input
           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) 
      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
           start_from_model=.false.
+          iranconf=(indpdb.le.0)
       endif
       if(start_from_model .and. (me.eq.king .or. .not. out1file))
      &    write(iout,*) 'START_FROM_MODELS is ON'
-      if(start_from_model .and. rest) then 
-        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
-         write(iout,*) 'START_FROM_MODELS is OFF'
-         write(iout,*) 'remove restart keyword from input'
-        endif
-      endif
+c      if(start_from_model .and. rest) then 
+c        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+c         write(iout,*) 'START_FROM_MODELS is OFF'
+c         write(iout,*) 'remove restart keyword from input'
+c        endif
+c      endif
+      if (start_from_model) nmodel_start=constr_homology
       if (homol_nset.gt.1)then
          call card_concat(controlcard)
          read(controlcard,*) (waga_homology(i),i=1,homol_nset) 
         tpl_k_rescore="template"//kic2//".sco"
 
         unres_pdb=.false.
+        nres_temp=nres
         if (read2sigma) then
           call readpdb_template(k)
         else
           call readpdb
         endif
+        nres_chomo(k)=nres
+        nres=nres_temp
 c
 c     Distance restraints
 c
 c          ... --> odl(k,ii)
 C Copy the coordinates from reference coordinates (?)
-        do i=1,2*nres
+        do i=1,2*nres_chomo(k)
           do j=1,3
             c(j,i)=cref(j,i)
 c           write (iout,*) "c(",j,i,") =",c(j,i)
@@ -3251,6 +3254,8 @@ c    &                       constr_homology
           enddo
         lim_odl=ii
         endif
+c        write (iout,*) "Distance restraints set"
+c        call flush(iout)
 c
 c     Theta, dihedral and SC retraints
 c
@@ -3291,6 +3296,8 @@ c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
           enddo
           lim_dih=nct-nnt-2 
         endif
+c        write (iout,*) "Dihedral angle restraints set"
+c        call flush(iout)
 
         if (waga_theta.gt.0.0d0) then
 c         open (ientin,file=tpl_k_sigma_theta,status='old')
@@ -3326,6 +3333,8 @@ c                             rescore(k,i-2) !  right expression ?
 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
           enddo
         endif
+c        write (iout,*) "Angle restraints set"
+c        call flush(iout)
 
         if (waga_d.gt.0.0d0) then
 c       open (ientin,file=tpl_k_sigma_d,status='old')
@@ -3361,51 +3370,66 @@ c              read (ientin,*) sigma_d(k,i) ! 1st variant
           enddo
         endif
       enddo
+c      write (iout,*) "SC restraints set"
+c      call flush(iout)
 c
 c remove distance restraints not used in any model from the list
 c shift data in all arrays
 c
+c      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
       if (waga_dist.ne.0.0d0) then
         ii=0
         liiflag=.true.
+        lfirst=.true.
         do i=nnt,nct-2 
          do j=i+2,nct 
           ii=ii+1
-          if (ii_in_use(ii).eq.0.and.liiflag) then
+c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+c     &            .and. distal.le.dist2_cut ) then
+c          write (iout,*) "i",i," j",j," ii",ii
+c          call flush(iout)
+          if (ii_in_use(ii).eq.0.and.liiflag.or.
+     &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
             liiflag=.false.
-            iistart=ii
+            i10=ii
+            if (lfirst) then
+              lfirst=.false.
+              iistart=ii
+            else
+              if(i10.eq.lim_odl) i10=i10+1
+              do ki=0,i10-i01-1
+               ires_homo(iistart+ki)=ires_homo(ki+i01)
+               jres_homo(iistart+ki)=jres_homo(ki+i01)
+               ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+               do k=1,constr_homology
+                odl(k,iistart+ki)=odl(k,ki+i01)
+                sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+                l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+               enddo
+              enddo
+              iistart=iistart+i10-i01
+            endif
           endif
-          if (ii_in_use(ii).ne.0.and..not.liiflag.or.
-     &                   .not.liiflag.and.ii.eq.lim_odl) then
-             if (ii.eq.lim_odl) then
-              iishift=ii-iistart+1
-             else
-              iishift=ii-iistart
-             endif
+          if (ii_in_use(ii).ne.0.and..not.liiflag) then
+             i01=ii
              liiflag=.true.
-             do ki=iistart,lim_odl-iishift
-              ires_homo(ki)=ires_homo(ki+iishift)
-              jres_homo(ki)=jres_homo(ki+iishift)
-              ii_in_use(ki)=ii_in_use(ki+iishift)
-              do k=1,constr_homology
-               odl(k,ki)=odl(k,ki+iishift)
-               sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
-               l_homo(k,ki)=l_homo(k,ki+iishift)
-              enddo
-             enddo
-             ii=ii-iishift
-             lim_odl=lim_odl-iishift
           endif
          enddo
         enddo
+        lim_odl=iistart-1
       endif
-
+c      write (iout,*) "Removing distances completed"
+c      call flush(iout)
       endif ! .not. klapaucjusz
 
       if (constr_homology.gt.0) call homology_partition
+c      write (iout,*) "After homology_partition"
+c      call flush(iout)
       if (constr_homology.gt.0) call init_int_table
-c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
-c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c      write (iout,*) "After init_int_table"
+c      call flush(iout)
+c      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 c
 c Print restraints
 c
@@ -3455,6 +3479,7 @@ c----------------------------------------------------------------------
 #endif
 c------------------------------------------------------------------------------
       subroutine copy_to_tmp(source)
+      implicit none
       include "DIMENSIONS"
       include "COMMON.IOUNITS"
       character*(*) source
@@ -3474,6 +3499,7 @@ c------------------------------------------------------------------------------
       end
 c------------------------------------------------------------------------------
       subroutine move_from_tmp(source)
+      implicit none
       include "DIMENSIONS"
       include "COMMON.IOUNITS"
       character*(*) source
@@ -3490,13 +3516,14 @@ c------------------------------------------------------------------------------
 C
 C Initialize random number generator
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
       logical OKRandom, prng_restart
       real*8  r1
       integer iseed_array(4)
+      integer error_msg,ierr
 #endif
       include 'COMMON.IOUNITS'
       include 'COMMON.TIME1'
@@ -3512,6 +3539,8 @@ C
       include 'COMMON.MD'
       include 'COMMON.FFIELD'
       include 'COMMON.SETUP'
+      integer i,iseed
+      double precision seed,ran_number
       iseed=-dint(dabs(seed))
       if (iseed.eq.0) then
         write (iout,'(/80(1h*)/20x,a/80(1h*))') 
@@ -3573,13 +3602,14 @@ c        r1 = prng_next(me)
       end
 c----------------------------------------------------------------------
       subroutine read_klapaucjusz
-
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
 #endif
       include 'COMMON.SETUP'
       include 'COMMON.CONTROL'
+      include 'COMMON.HOMOLOGY'
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
       include 'COMMON.MD'
@@ -3588,13 +3618,16 @@ c----------------------------------------------------------------------
       include 'COMMON.NAMES'
       character*256 fragfile
       integer ninclust(maxclust),inclust(max_template,maxclust),
-     &  nresclust(maxclust),iresclust(maxres,maxclust)
+     &  nresclust(maxclust),iresclust(maxres,maxclust),nclust
 
       character*2 kic2
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
-      integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
+      integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
+     & ik,ll,ii,kk,iistart,iishift,lim_xx
+      double precision distal
       logical lprn /.true./
+      integer nres_temp
       integer ilen
       external ilen
       logical liiflag
@@ -3629,7 +3662,10 @@ c Read pdb files
         stop
   34    continue
         unres_pdb=.false.
+        nres_temp=nres
         call readpdb_template(k)
+        nres_chomo(k)=nres
+        nres=nres_temp
         do i=1,nres
           rescore(k,i)=0.2d0
           rescore2(k,i)=1.0d0
@@ -3663,6 +3699,8 @@ c     Distance restraints
 c
 c          ... --> odl(k,ii)
 C Copy the coordinates from reference coordinates (?)
+        nres_temp=nres
+        nres=nres_chomo(k)
         do i=1,2*nres
           do j=1,3
             c(j,i)=chomo(j,i,k)
@@ -3670,11 +3708,12 @@ c           write (iout,*) "c(",j,i,") =",c(j,i)
           enddo
         enddo
         call int_from_cart(.true.,.false.)
-        call sc_loc_geom(.true.)
+        call sc_loc_geom(.false.)
         do i=1,nres
           thetaref(i)=theta(i)
           phiref(i)=phi(i)
         enddo
+        nres=nres_temp
         if (waga_dist.ne.0.0d0) then
           ii=0
           do i = nnt,nct-2