make cp src-HCD-5D
[unres.git] / source / unres / src-HCD-5D / readrtns_CSA.F
index 66f7c17..4765a41 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,25 +152,24 @@ 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 reada(controlcard,"T_BATH",t_bath,300.0d0)
       minim=(index(controlcard,'MINIMIZE').gt.0)
@@ -297,9 +304,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 +342,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 +355,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 +363,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 +378,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 +444,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 +468,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)
@@ -545,7 +551,7 @@ c  if performing umbrella sampling, fragments constrained are read from the frag
      &    "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 +687,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 +704,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 +720,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
+      double precision sumv
 C
 C Read PDB structure if applicable
 C
@@ -789,25 +801,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)
@@ -840,10 +833,18 @@ c      print '(20i4)',(itype(i),i=1,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
+      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,8 +856,11 @@ 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
@@ -1372,10 +1376,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,6 +1397,7 @@ 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
@@ -1447,8 +1453,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 +1466,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 +1475,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 +1506,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 +1522,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 +1553,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 +1569,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 +1589,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'
@@ -1602,8 +1610,9 @@ C Generate CA distance constraints.
       include 'COMMON.DBASE'
       include 'COMMON.THREAD'
       include 'COMMON.TIME1'
-      dimension itype_pdb(maxres)
+      integer i,j,itype_pdb(maxres)
       common /pizda/ itype_pdb
+      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,
@@ -1638,10 +1647,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 +1694,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 +1756,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 +1824,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 +1850,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 +1951,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 +1963,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,13 +2352,21 @@ 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
@@ -2461,9 +2382,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 +2393,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 +2445,7 @@ c-------------------------------------------------------------------------------
       end
 C---------------------------------------------------------------------------
       subroutine read_afminp
-            implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include 'mpif.h'
@@ -2533,7 +2456,8 @@ 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)
@@ -2546,22 +2470,24 @@ CCCC NOW PROPERTIES FOR AFM
         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 +2544,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 +2555,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 +2945,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,7 +2973,8 @@ 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,iishift
       integer ilen
       external ilen
       logical liiflag
@@ -3055,6 +2986,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
@@ -3455,6 +3387,7 @@ c----------------------------------------------------------------------
 #endif
 c------------------------------------------------------------------------------
       subroutine copy_to_tmp(source)
+      implicit none
       include "DIMENSIONS"
       include "COMMON.IOUNITS"
       character*(*) source
@@ -3474,6 +3407,7 @@ c------------------------------------------------------------------------------
       end
 c------------------------------------------------------------------------------
       subroutine move_from_tmp(source)
+      implicit none
       include "DIMENSIONS"
       include "COMMON.IOUNITS"
       character*(*) source
@@ -3490,13 +3424,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 +3447,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 +3510,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,12 +3526,14 @@ 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 ilen
       external ilen
@@ -3670,7 +3610,7 @@ 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)