update
[unres.git] / source / unres / src_MD-M / gen_rand_conf.F
index d870f55..6bb4f1a 100644 (file)
@@ -1,4 +1,4 @@
-      subroutine gen_rand_conf(nstart,*)
+      subroutine gen_rand_conf(nstart,nend,*)
 C Generate random conformation or chain cut and regrowth.
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
@@ -10,15 +10,21 @@ C Generate random conformation or chain cut and regrowth.
       include 'COMMON.MCM'
       include 'COMMON.GEO'
       include 'COMMON.CONTROL'
-      logical overlap,back,fail
-cd    print *,' CG Processor',me,' maxgen=',maxgen
+      logical overlap,back,fail,ldir
+      maxgen=10000
+      write (iout,*) ' CG Processor',me,' maxgen=',maxgen
       maxsi=100
-cd    write (iout,*) 'Gen_Rand_conf: nstart=',nstart
+      ldir = nstart.le.nend
+      write (iout,*) 'Gen_Rand_conf: nstart=',nstart,' nend',nend,
+     & ' ldir',ldir
+
+      if (ldir) then
+
       if (nstart.lt.5) then
-        it1=itype(2)
-        phi(4)=gen_phi(4,itype(2),itype(3))
+        it1=iabs(itype(2))
+        phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3)))
 c       write(iout,*)'phi(4)=',rad2deg*phi(4)
-        if (nstart.lt.3) theta(3)=gen_theta(itype(2),pi,phi(4))
+        if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4))
 c       write(iout,*)'theta(3)=',rad2deg*theta(3) 
         if (it1.ne.10) then
           nsi=0
@@ -42,7 +48,7 @@ c       write(iout,*)'theta(3)=',rad2deg*theta(3)
       nit=0
       niter=0
       back=.false.
-      do while (i.le.nres .and. niter.lt.maxgen)
+      do while (i.le.nend .and. niter.lt.maxgen)
         if (i.lt.nstart) then
           if(iprint.gt.1) then
           write (iout,'(/80(1h*)/2a/80(1h*))') 
@@ -54,9 +60,9 @@ c       write(iout,*)'theta(3)=',rad2deg*theta(3)
           endif
           return1
         endif
-       it1=itype(i-1)
-       it2=itype(i-2)
-       it=itype(i)
+       it1=iabs(itype(i-1))
+       it2=iabs(itype(i-2))
+       it=iabs(itype(i))
 c       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
 c    &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
        phi(i+1)=gen_phi(i+1,it1,it)
@@ -64,7 +70,7 @@ c    &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
           phi(i)=gen_phi(i+1,it2,it1)
 c         print *,'phi(',i,')=',phi(i)
          theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
-         if (it2.ne.10) then
+         if (it2.ne.10 .and. it2.ne.ntyp1) then
             nsi=0
             fail=.true.
             do while (fail.and.nsi.le.maxsi)
@@ -76,7 +82,7 @@ c         print *,'phi(',i,')=',phi(i)
          call locate_next_res(i-1)
        endif
        theta(i)=gen_theta(it1,phi(i),phi(i+1))
-        if (it1.ne.10) then 
+        if (it1.ne.10 .and. it1.ne.ntyp1) then 
         nsi=0
         fail=.true.
         do while (fail.and.nsi.le.maxsi)
@@ -86,7 +92,7 @@ c         print *,'phi(',i,')=',phi(i)
         if (nsi.gt.maxsi) return1
         endif
        call locate_next_res(i)
-       if (overlap(i-1)) then
+       if (overlap(i-1,ldir)) then
          if (nit.lt.maxnit) then
            back=.true.
            nit=nit+1
@@ -121,23 +127,134 @@ c         print *,'phi(',i,')=',phi(i)
        c(j,nres+1)=c(j,1)
        c(j,nres+nres)=c(j,nres)
       enddo
+
+      else
+
+      maxnit=100
+
+      i=nstart
+      nit=0
+      niter=0
+      back=.false.
+      do while (i.ge.nend .and. niter.lt.maxgen)
+        if (i.gt.nstart) then
+          if(iprint.gt.1) then
+          write (iout,'(/80(1h*)/2a/80(1h*))') 
+     &          'Generation procedure went down to ',
+     &          'chain beginning. Cannot continue...'
+          write (*,'(/80(1h*)/2a/80(1h*))') 
+     &          'Generation procedure went down to ',
+     &          'chain beginning. Cannot continue...'
+          endif
+          return1
+        endif
+       it1=iabs(itype(i+1))
+       it2=iabs(itype(i+2))
+       it3=iabs(itype(i+3))
+c       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
+c    &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
+       phi(i+2)=gen_phi(i+2,it1,it2)
+       if (back) then
+          phi(i+3)=gen_phi(i+3,it2,it3)
+          phi(i+4)=gen_phi(i+3,it2,it3)
+c         print *,'phi(',i+3,')=',phi(i+3)
+         theta(i+3)=gen_theta(it2,phi(i+3),phi(i+4))
+         if (it2.ne.10 .and. it2.ne.ntyp1) then
+            nsi=0
+            fail=.true.
+            do while (fail.and.nsi.le.maxsi)
+              call gen_side(it2,theta(i+3),alph(i+2),omeg(i+2),fail)
+              nsi=nsi+1
+            enddo
+            if (nsi.gt.maxsi) return1
+          endif
+         call locate_prev_res(i+1)
+       endif
+       theta(i+2)=gen_theta(it1,phi(i+2),phi(i+3))
+        write (iout,*) "i"," theta",theta(i+2)," phi",phi(i+2),phi(i+3)
+        if (it1.ne.10 .and. it1.ne.ntyp1) then 
+        nsi=0
+        fail=.true.
+        do while (fail.and.nsi.le.maxsi)
+          call gen_side(it1,theta(i+2),alph(i+1),omeg(i+1),fail)
+          nsi=nsi+1
+        enddo
+        write (iout,*) "i",i," alpha",alph(i+1)," omeg",omeg(i+1),
+     &     " fail",fail
+        if (nsi.gt.maxsi) return1
+        endif
+       call locate_prev_res(i)
+        write (iout,*) "After locate_prev_res i=",i
+        write (iout,'(3f10.5,5x,3f10.5)') (c(l,i),l=1,3),
+     &    (c(l,i+nres),l=1,3)
+        write (iout,'(3f10.5,5x,3f10.5)') (c(l,i+1),l=1,3),
+     &    (c(l,i+1+nres),l=1,3)
+        call flush(iout)
+       if (overlap(i+1,ldir)) then
+          write (iout,*) "OVERLAP",nit
+          call flush(iout)
+         if (nit.lt.maxnit) then
+           back=.true.
+           nit=nit+1
+          else
+           nit=0
+           if (i.gt.3) then
+             back=.true.
+             i=i+1
+            else
+             write (iout,'(a)') 
+     &  'Cannot generate non-overlaping conformation. Increase MAXNIT.'
+             write (*,'(a)') 
+     &  'Cannot generate non-overlaping conformation. Increase MAXNIT.'
+             return1
+           endif
+          endif
+        else
+          write (iout,*) "================>> NO OVERLAP",i
+          call flush(iout)
+         back=.false.
+         nit=0 
+         i=i-1
+        endif
+       niter=niter+1
+        write (iout,*) "iter",niter
+      enddo
+      if (niter.ge.maxgen) then
+       write (iout,'(a,2i5)') 
+     &'Too many trials in backward conformation generation',niter,maxgen
+       write (*,'(a,2i5)') 
+     &'Too many trials in backward conformation generation',niter,maxgen
+       return1
+      endif
+
+      endif
+
+      write (iout,*) "HERE!!!!"
+
       return
       end
 c-------------------------------------------------------------------------
-      logical function overlap(i)
+      logical function overlap(i,ldir)
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
       include 'COMMON.CHAIN'
       include 'COMMON.INTERACT'
       include 'COMMON.FFIELD'
       data redfac /0.5D0/
+      logical ldir
+      logical lprn /.false./
       overlap=.false.
-      iti=itype(i)
+      iti=iabs(itype(i))
       if (iti.gt.ntyp) return
 C Check for SC-SC overlaps.
 cd    print *,'nnt=',nnt,' nct=',nct
+c      write (iout,*) "overlap i",i," ldir",ldir
+
+      if (ldir) then
+
       do j=nnt,i-1
-        itj=itype(j)
+        itj=iabs(itype(j))
         if (j.lt.i-1 .or. ipot.ne.4) then
           rcomp=sigmaii(iti,itj)
         else 
@@ -146,9 +263,9 @@ cd    print *,'nnt=',nnt,' nct=',nct
 cd      print *,'j=',j
        if (dist(nres+i,nres+j).lt.redfac*rcomp) then
           overlap=.true.
-c        print *,'overlap, SC-SC: i=',i,' j=',j,
-c     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
-c     &     rcomp
+          if (lprn) write (iout,*) 'overlap, SC-SC: i=',i,' j=',j,
+     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
+     &     rcomp
          return
         endif
       enddo
@@ -159,10 +276,10 @@ C SCs.
        c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
       enddo
       do j=nnt,i-2
-       itj=itype(j)
-cd      print *,'overlap, p-Sc: i=',i,' j=',j,
-cd   &         ' dist=',dist(nres+j,maxres2+1)
+       itj=iabs(itype(j))
        if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
+          if (lprn) write (iout,*) 'overlap, p-Sc: i=',i,' j=',j,
+     &         ' dist=',dist(nres+j,maxres2+1)
          overlap=.true.
          return
         endif
@@ -173,9 +290,9 @@ C groups.
        do k=1,3
          c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1))
         enddo
-cd      print *,'overlap, SC-p: i=',i,' j=',j,
-cd   &         ' dist=',dist(nres+i,maxres2+1)
        if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then
+          if (lprn) write (iout,*) 'overlap, SC-p: i=',i,' j=',j,
+     &         ' dist=',dist(nres+i,maxres2+1)
           overlap=.true.
          return
         endif
@@ -189,26 +306,104 @@ C Check for p-p overlaps
        do k=1,3
          c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
         enddo
-cd      print *,'overlap, p-p: i=',i,' j=',j,
-cd   &         ' dist=',dist(maxres2+1,maxres2+2)
         if(iteli.ne.0.and.itelj.ne.0)then
         if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then
+          if (lprn) write (iout,*) 'overlap, p-p: i=',i,' j=',j,
+     &         ' dist=',dist(maxres2+1,maxres2+2)
+          overlap=.true.
+          return
+        endif
+        endif
+      enddo
+
+      else
+
+      if (lprn) write (iout,*) "start",i+1," end",nres_start+nsup
+      do j=i+1,nres_start+nsup
+        itj=iabs(itype(j))
+        if (j.gt.i+1 .or. ipot.ne.4) then
+          rcomp=sigmaii(iti,itj)
+        else 
+          rcomp=sigma(iti,itj)
+        endif
+cd      print *,'j=',j
+c        write (iout,*) "SCSC j",j," dist",dist(nres+i,nres+j),
+c    &   " radfac",redfac*rcomp
+       if (dist(nres+i,nres+j).lt.redfac*rcomp) then
+          overlap=.true.
+          if (lprn) write (iout,*) 'overlap, SC-SC: i=',i,' j=',j,
+     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
+     &     rcomp
+         return
+        endif
+      enddo
+C Check for overlaps between the added peptide group and the succeeding
+C SCs.
+      iteli=itel(i)
+      do j=1,3
+       c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+      enddo
+      do j=i+2,nstart_res+nsup
+       itj=iabs(itype(j))
+cd        write (iout,*) "pSC j",j," dist",dist(nres+j,maxres2+1)
+       if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
+          if (lprn) write (iout,*) 'overlap, p-Sc: i=',i,' j=',j,
+     &         ' dist=',dist(nres+j,maxres2+1)
+         overlap=.true.
+         return
+        endif
+      enddo
+C Check for overlaps between the added side chain and the succeeding peptide
+C groups.
+      do j=i+2,nstart_seq+nsup
+       do k=1,3
+         c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1))
+        enddo
+cd        write (iout,*) "SCp j",j," dist",dist(nres+i,maxres2+1)
+       if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then
+          if (lprn) write (iout,*) 'overlap, SC-p: i=',i,' j=',j,
+     &         ' dist=',dist(nres+i,maxres2+1)
+          overlap=.true.
+         return
+        endif
+      enddo
+C Check for p-p overlaps
+      do j=1,3
+       c(j,maxres2+2)=0.5D0*(c(j,i)+c(j,i+1))
+      enddo
+      do j=i+2,nres_start+nsup
+        itelj=itel(j)
+       do k=1,3
+         c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
+        enddo
+        if(iteli.ne.0.and.itelj.ne.0)then
+cd        write (iout,*) "pp j",j," dist",dist(maxres2+1,maxres2+2)
+        if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then
+           if (lprn) write (iout,*) 'overlap, p-p: i=',i,' j=',j,
+     &         ' dist=',dist(maxres2+1,maxres2+2)
           overlap=.true.
           return
         endif
         endif
       enddo
+
+      endif 
+
       return
       end
 c--------------------------------------------------------------------------
       double precision function gen_phi(i,it1,it2)
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
+      include "COMMON.TORCNSTR"
       include 'COMMON.GEO'
       include 'COMMON.BOUNDS'
-c      gen_phi=ran_number(-pi,pi) 
+      if (raw_psipred .or. ndih_constr.eq.0) then
+        gen_phi=ran_number(-pi,pi) 
+      else
 C 8/13/98 Generate phi using pre-defined boundaries
-      gen_phi=ran_number(phibound(1,i),phibound(2,i)) 
+        gen_phi=ran_number(phibound(1,i),phibound(2,i)) 
+      endif
       return
       end
 c---------------------------------------------------------------------------
@@ -238,7 +433,8 @@ c     print *,'gen_theta: it=',it
       endif  
       thet_pred_mean=a0thet(it)
       do k=1,2
-        thet_pred_mean=thet_pred_mean+athet(k,it)*y(k)+bthet(k,it)*z(k)
+        thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k)
+     &     +bthet(k,it,1,1)*z(k)
       enddo
       sig=polthet(3,it)
       do j=2,0,-1
@@ -779,7 +975,7 @@ c     overlapping residues left, or false otherwise (success)
 
         do ires=1,ioverlap_last 
           i=ioverlap(ires)
-          iti=itype(i)
+          iti=iabs(itype(i))
           if (iti.ne.10) then
             nsi=0
             fail=.true.
@@ -791,7 +987,7 @@ c     overlapping residues left, or false otherwise (success)
           endif
         enddo
 
-        call chainbuild
+        call chainbuild_extconf
         call overlap_sc_list(ioverlap,ioverlap_last)
 c        write (iout,*) 'Overlaping residues ',ioverlap_last,
 c     &           (ioverlap(j),j=1,ioverlap_last)
@@ -839,8 +1035,8 @@ C Check for SC-SC overlaps and mark residues
 c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -852,7 +1048,7 @@ c
        do iint=1,nint_gr(i)
          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
-            itypj=itype(j)
+            itypj=iabs(itype(j))
             dscj_inv=dsc_inv(itypj)
             sig0ij=sigma(itypi,itypj)
             chi1=chi(itypi,itypj)