corrections
[unres.git] / source / unres / src-HCD-5D / readrtns_CSA.F
index 16c0f37..d76b29e 100644 (file)
@@ -187,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"
@@ -296,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
@@ -736,7 +741,7 @@ C
       integer ilen
       external ilen
       integer iperm,tperm
-      integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2
+      integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
       double precision sumv
 C
 C Read PDB structure if applicable
@@ -831,7 +836,7 @@ 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
@@ -844,7 +849,7 @@ cd      print *,'NNT=',NNT,' NCT=',NCT
         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
+      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
@@ -870,9 +875,9 @@ c      enddo
       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
@@ -1097,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
@@ -1163,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
 
 
@@ -1172,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)
@@ -1405,7 +1454,13 @@ C Read information about disulfide bridges.
       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.
@@ -1614,9 +1669,11 @@ C Generate CA distance constraints.
       include 'COMMON.CONTROL'
       include 'COMMON.DBASE'
       include 'COMMON.THREAD'
+      include 'COMMON.SPLITELE'
       include 'COMMON.TIME1'
       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
@@ -1627,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:' 
@@ -2375,10 +2435,10 @@ c------------------------------------------------------------------------------
       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
@@ -2468,7 +2528,7 @@ c      print *, "wchodze"
       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
@@ -2979,7 +3039,7 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
       integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
-     & ik,iistart
+     & ik,iistart,nres_temp
       integer ilen
       external ilen
       logical liiflag,lfirst
@@ -3016,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)
@@ -3562,6 +3627,7 @@ c----------------------------------------------------------------------
      & ik,ll,ii,kk,iistart,iishift,lim_xx
       double precision distal
       logical lprn /.true./
+      integer nres_temp
       integer ilen
       external ilen
       logical liiflag
@@ -3596,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
@@ -3630,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)
@@ -3642,6 +3713,7 @@ c           write (iout,*) "c(",j,i,") =",c(j,i)
           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