Adam's changes to wham and cluster following previous commit
[unres.git] / source / cluster / wham / src / readrtns.F
index e68a9d3..037eb1f 100644 (file)
@@ -37,7 +37,8 @@ C
       min_var=(index(controlcard,'MINVAR').gt.0)
       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
-      call readi(controlcard,'NCUT',ncut,1)
+      call readi(controlcard,'NCUT',ncut,0)
+      call readi(controlcard,'NCLUST',nclust,5)
       call readi(controlcard,'NSTART',nstart,0)
       call readi(controlcard,'NEND',nend,0)
       call reada(controlcard,'ECUT',ecut,10.0d0)
@@ -46,7 +47,8 @@ C
       lgrp=(index(controlcard,'LGRP').gt.0)
       caonly=(index(controlcard,'CA_ONLY').gt.0)
       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
-      call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
+      if (ncut.gt.0) 
+     & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
       call readi(controlcard,'IOPT',iopt,2) 
       lside = index(controlcard,"SIDE").gt.0
       efree = index(controlcard,"EFREE").gt.0
@@ -69,8 +71,17 @@ C
       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
       write (iout,*) "with_homology_constr ",with_dihed_constr,
      & " CONSTR_HOMOLOGY",constr_homology
+      print_homology_restraints=
+     & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
+      print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
+      print_homology_models=
+     & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
 
+#ifdef AIX
+      call flush_(iout)
+#else
       call flush(iout)
+#endif
       if (min_var) iopt=1
       return
       end
@@ -247,7 +258,11 @@ C 12/1/95 Added weight for the multi-body term WCORR
         bad(i,2)=scalscp*bad(i,2)
       enddo
 
+#ifdef AIX
+      call flush_(iout)
+#else
       call flush(iout)
+#endif
       print *,'indpdb=',indpdb,' pdbref=',pdbref
 
 C Read sequence if not taken from the pdb file.
@@ -260,6 +275,19 @@ C Convert sequence to numeric code
       do i=1,nres
         itype(i)=rescode(i,sequence(i),iscode)
       enddo
+      if (itype(2).eq.10) then
+        write (iout,*)
+     &   "Glycine is the first full residue, initial dummy deleted"
+        do i=1,nres
+          itype(i)=itype(i+1)
+        enddo
+        nres=nres-1
+      endif
+      if (itype(nres).eq.10) then
+        write (iout,*)
+     &   "Glycine is the last full residue, terminal dummy deleted"
+        nres=nres-1
+      endif
       print *,nres
       print '(20i4)',(itype(i),i=1,nres)
 
@@ -325,18 +353,36 @@ C     both routines should be in dfa.F file!!
 
       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
+#ifdef DEBUG
        write (iout,*) "Calling init_dfa_vars"
-       call flush(iout)
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
+#endif
        call init_dfa_vars
+#ifdef DEBUG
        write (iout,*) 'init_dfa_vars finished!'
-       call flush(iout)
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
+#endif
        call read_dfa_info
+#ifdef DEBUG
        write (iout,*) 'read_dfa_info finished!'
-       call flush(iout)
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
+#endif
       endif
 
       if (constr_homology.gt.0) then
-        call read_constr_homology
+        call read_constr_homology(print_homology_restraints)
       endif
 
 c      if (pdbref) then
@@ -406,7 +452,7 @@ c      endif
             enddo
           enddo
         endif
-        call contact(.true.,ncont_ref,icont_ref)
+        call contact(print_contact_map,ncont_ref,icont_ref)
       endif
 c Read distance restraints
       if (constr_dist.gt.0) then
@@ -751,7 +797,11 @@ c      write (iout,'(a)') controlcard
       do i=1,npair_
         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
       enddo
+#ifdef AIX
+      call flush_(iout)
+#else
       call flush(iout)
+#endif
       if (.not.refstr .and. nfrag_.gt.0) then
         write (iout,*) 
      &  "ERROR: no reference structure to compute distance restraints"
@@ -764,13 +814,17 @@ c      write (iout,'(a)') controlcard
      &   " but distance restraints between pairs requested"
         stop 
       endif 
+#ifdef AIX
+      call flush_(iout)
+#else
       call flush(iout)
+#endif
       do i=1,nfrag_
         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
      &    ifrag_(2,i)=nstart_sup+nsup-1
 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
-        call flush(iout)
+c        call flush(iout)
         if (wfrag_(i).gt.0.0d0) then
         do j=ifrag_(1,i),ifrag_(2,i)-1
           do k=j+1,ifrag_(2,i)
@@ -842,12 +896,16 @@ c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
       enddo
+#ifdef AIX
+      call flush_(iout)
+#else
       call flush(iout)
+#endif
       return
       end
 
 c====-------------------------------------------------------------------
-      subroutine read_constr_homology
+      subroutine read_constr_homology(lprn)
 
       include 'DIMENSIONS'
 #ifdef MPI
@@ -875,7 +933,8 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
       integer ki, i, j, k, l
-      logical lprn /.true./
+      logical lprn
+      logical unres_pdb
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
 c
@@ -890,8 +949,14 @@ c -----------------------------------------------------------------
 c
 c
 c Alternative: reading from input
+#ifdef DEBUG
       write (iout,*) "BEGIN READ HOMOLOGY INFO"
+#ifdef AIX
+      call flush_(iout)
+#else
       call flush(iout)
+#endif
+#endif
       call card_concat(controlcard)
       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
@@ -901,7 +966,7 @@ c Alternative: reading from input
 
       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
       if (homol_nset.gt.1)then
-         call readi(controlcard,"ISET",iset,homol_nset)
+         call readi(controlcard,"ISET",iset,1)
          call card_concat(controlcard)
          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
       else
@@ -909,11 +974,15 @@ c Alternative: reading from input
         waga_homology(1)=1.0
       endif
 c
+#ifdef DEBUG
       write(iout,*) "read_constr_homology iset",iset
       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
+#ifdef AIX
+      call flush_(iout)
+#else
       call flush(iout)
-
-
+#endif
+#endif
 cd      write (iout,*) "nnt",nnt," nct",nct
 cd      call flush(iout)
 
@@ -942,7 +1011,7 @@ c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
       do k=1,constr_homology
 
         read(inp,'(a)') pdbfile
-        write (iout,*) "k ",k," pdbfile ",pdbfile
+c        write (iout,*) "k ",k," pdbfile ",pdbfile
 c  Next stament causes error upon compilation (?)
 c       if(me.eq.king.or. .not. out1file)
 c         write (iout,'(2a)') 'PDB data will be read from file ',
@@ -977,7 +1046,6 @@ c       tpl_k_sigma_d="template"//kic2//".sigma_d"
           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
      &      (crefjlee(j,i+nres),j=1,3)
         enddo
-#endif
         write (iout,*) "READ HOMOLOGY INFO"
         write (iout,*) "read_constr_homology x: after reading pdb file"
         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
@@ -986,7 +1054,12 @@ c       tpl_k_sigma_d="template"//kic2//".sigma_d"
         write (iout,*) "waga_theta",waga_theta
         write (iout,*) "waga_d",waga_d
         write (iout,*) "dist_cut",dist_cut
-        call flush(iout)
+#endif
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
 
 c
 c     Distance restraints