X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fgen_rand_conf_mchain.F;fp=source%2Funres%2Fsrc-HCD-5D%2Fgen_rand_conf_mchain.F;h=4614e8744b5c6eb8b8f2b306973a1830ca949292;hb=43dfdc33f5b338546258b5a882326f97a76b5bd4;hp=0000000000000000000000000000000000000000;hpb=68aa1d135805a14583695611acfcf637748b4466;p=unres.git

diff --git a/source/unres/src-HCD-5D/gen_rand_conf_mchain.F b/source/unres/src-HCD-5D/gen_rand_conf_mchain.F
new file mode 100644
index 0000000..4614e87
--- /dev/null
+++ b/source/unres/src-HCD-5D/gen_rand_conf_mchain.F
@@ -0,0 +1,424 @@
+      subroutine gen_rand_conf_mchain(nstart0,*)
+C Generate random conformation or chain cut and regrowth.
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.MCM'
+      include 'COMMON.GEO'
+      include 'COMMON.CONTROL'
+      logical overlap_mchain,back,fail
+      integer nstart,nstart0
+      integer i,ii,iii,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit
+      integer igr,iequi,ichain,nnres
+      double precision aux
+      double precision gen_theta,gen_phi,dist,ran_number,scalar
+c      write (iout,*)  'gen_rand_conf_mchain: maxgen=',maxgen
+      nstart=nstart0
+      maxsi=100
+c      write (iout,*) 'Gen_Rand_conf_mchain: nstart=',nstart,
+c     &   " nchain_group",nchain_group
+
+      DO IGR=1,NCHAIN_GROUP
+
+      DO IEQUI=1,NEQUIV(IGR)
+
+      ichain=iequiv(iequi,igr)
+
+      i=chain_border1(1,ichain)+nstart-1
+      if (nstart.eq.1) then
+        do j=1,3
+          c(j,i)=ran_number(-15.0d0,15.0d0)
+          dc(j,i-1)=c(j,i)
+        enddo
+      endif
+      if (nstart.le.2) then
+
+      do j=1,3
+        dc_norm(j,i)=ran_number(-1.0d0,1.0d0)
+      enddo
+      aux=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i)))
+      do j=1,3
+        dc_norm(j,i)=dc_norm(j,i)/aux
+      enddo
+      if (itype(i).eq.ntyp1) then
+        do j=1,3
+          dc(j,i)=1.9d0*dc_norm(j,i)
+        enddo
+      else 
+        do j=1,3
+          dc(j,i)=3.8d0*dc_norm(j,i)
+        enddo
+      endif
+      do j=1,3
+        c(j,i+1)=c(j,i)+dc(j,i)
+      enddo
+      endif
+c      if (nstart.lt.5) then
+      if (nstart.le.2) then
+      it1=iabs(itype(i+2))
+      phi(i+3)=gen_phi(i+3,iabs(itype(i+1)),iabs(itype(i+2)))
+c       write(iout,*)'phi(4)=',rad2deg*phi(4)
+      theta(i+2)=gen_theta(iabs(itype(i+2)),pi,phi(i+3))
+      if (it1.ne.10) 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
+        if (nsi.gt.maxsi) then
+          write (iout,'(a,i7,a,i7,a,i7,a,i7)') 
+     &    'Problem in SC rotamer generation, residue',ii,
+     &    ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
+          return1
+        endif
+      endif ! it1.ne.10
+      call orig_frame_chain(i)
+      else
+c      call orig_frame_chain(i-1)
+c      write (iout,*) "calling refsys",i
+      call refsys(i,i-1,i-2,prod(1,1,i-1),
+     &  prod(1,2,i-1),prod(1,3,i-1),fail)
+c      write (iout,*) "after refsys",i
+#ifdef DEBUG
+      write (iout,*) "dc_norm(:",i-1,") and prod"
+      do j=1,3
+        write (iout,*) j,dc_norm(j,i-1),(prod(j,k,i-1),k=1,3)
+      enddo
+#endif
+      endif
+
+      ENDDO
+
+      nstart=nstart+1
+      ii=nstart
+
+      maxnit=5000
+
+      nit=0
+      niter=0
+      back=.false.
+      nnres=chain_border1(2,iequiv(1,igr))-
+     &  chain_border1(1,iequiv(1,igr))+1
+#ifdef DEBUG
+      write (iout,*) "chain group",igr," chains",
+     &    (iequiv(j,igr),j=1,nequiv(igr))
+      write (iout,*) "ii",ii," nnres",nnres
+#endif
+      do while (ii.le. nnres .and. niter.lt.maxgen)
+
+        ichain=iequiv(1,igr)
+        i=ii-1+chain_border1(1,ichain)
+#ifdef DEBUG
+        write (iout,*) "ii",ii," nnres",nnres," ichain",ichain," i",i,
+     &    "niter",niter," back",back," nstart",nstart
+#endif
+        if (ii.lt.nstart) then
+c          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...'
+c          endif
+          return1
+        endif
+        it1=iabs(itype(i-1))
+        it2=iabs(itype(i-2))
+        it=iabs(itype(i))
+        if (it.eq.ntyp1 .and. it1.eq.ntyp1) then
+          vbld(i)=ran_number(3.8d0,10.0d0)
+          vbld_inv(i)=1.0d0/vbld(i)
+        endif
+#ifdef DEBUG
+        write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,
+     &    ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen
+#endif
+        phi(i)=gen_phi(i,it2,it1)
+        phi(i+1)=gen_phi(i+1,it1,it)
+#ifdef DEBUG
+        write (iout,*) "phi",i,phi(i)," phi",i+1,phi(i+1)
+#endif
+        do iequi=2,nequiv(igr)
+          iii=ii+chain_border1(1,iequiv(iequi,igr))
+          phi(iii)=phi(i+1)
+        enddo
+#ifdef CHUJ
+        if (back) then
+          phi(i)=gen_phi(i+1,it2,it1)
+#ifdef DEBUG
+          write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it
+#endif
+          theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
+          if (theta(i-1).gt.2.68780478d0) theta(i-1)=2.68780478d0
+          do iequi=2,nequiv(igr)
+            iii=ii-1+chain_border1(1,iequiv(iequi,igr))
+            phi(iii)=phi(i+1)
+            theta(iii-1)=theta(i-1)
+          enddo
+          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-1),alph(i-2),omeg(i-2),fail)
+              nsi=nsi+1
+            enddo
+            do iequi=2,nequiv(igr)
+              iii=ii-3+chain_border1(1,iequiv(iequi,igr))
+              alph(iii)=alph(i-2)
+              omeg(iii)=omeg(i-2)
+            enddo
+#ifdef DEBUG
+            write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
+#endif
+            if (nsi.gt.maxsi) then
+              write (iout,'(a,i7,a,i7,a,i7,a,i7)') 
+     &        'Problem in SC rotamer generation, residue',ii,
+     &        ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
+              return1
+            endif
+
+          endif
+c          call locate_next_res(i-1)
+        endif
+#endif
+        if (it1.ne.ntyp1) then
+          theta(i)=gen_theta(it1,phi(i),phi(i+1))
+          if (theta(i).gt.2.68780478d0) theta(i)=2.68780478d0
+        else
+          theta(i)=ran_number(1.326d0,2.548d0)
+        endif
+#ifdef DEBUG
+        write (iout,*) "ii",ii," i",i," it1",it1,
+     &   " theta",theta(i)," phi",
+     &    phi(i),phi(i+1)
+#endif
+        if (it1.ne.10 .and. it1.ne.ntyp1) then 
+        nsi=0
+        fail=.true.
+        do while (fail.and.nsi.le.maxsi)
+c          theta(i)=gen_theta(it1,phi(i),phi(i+1))
+          call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
+          nsi=nsi+1
+        enddo
+#ifdef DEBUG
+        write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
+#endif
+c        write (iout,*) "After forward SC generation:",nsi,maxsi
+        if (nsi.gt.maxsi) then
+          write (iout,'(a,i7,a,i7,a,i7,a,i7)') 
+     &    'Problem in SC rotamer generation, residue',ii,
+     &    ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
+          return1
+        endif
+
+        endif
+
+        do iequi=2,nequiv(igr)
+          iii=ii-1+chain_border1(1,iequiv(iequi,igr))
+          theta(iii)=theta(i)
+          phi(iii)=phi(i)
+          alph(iii-1)=alph(i-1)
+          omeg(iii-1)=omeg(i-1)
+        enddo
+
+        DO IEQUI=1,NEQUIV(IGR)
+
+        ichain=iequiv(iequi,igr)
+        i=ii-1+chain_border1(1,ichain)
+#ifdef CHUJ
+        if (back) call locate_next_res(i-1)
+#endif
+        call locate_next_res(i)
+#ifdef DEBUG
+        write (iout,*) "i",i," ii",ii," ichain",ichain
+        write (iout,*) theta(i)*rad2deg,phi(i)*rad2deg,
+     &    alph(i-1)*rad2deg,omeg(i-1)*rad2deg
+        write (iout,*) (c(j,i),j=1,3),(c(j,i+nres-1),j=1,3)
+#endif
+        if (overlap_mchain(i-1,ii-1,ichain,igr)) then
+#ifdef DEBUG
+          write (iout,*) "***********overlap",i-1," nit",nit
+#endif
+          if (nit.lt.maxnit) then
+            back=.true.
+            nit=nit+1
+            exit
+          else
+#ifdef DEBUG
+            write (iout,*) "***********overlap maxnit exceeded",nit
+#endif
+            nit=0
+            if (ii.gt.3) then
+              back=.true.
+              ii=ii-1
+c              write (iout,*) "ii",ii
+              exit
+            else
+              write (iout,'(a,i7,a,i7,a,i7,a,i7)') 
+     &        'Cannot generate non-overlaping conformation, residue',ii,
+     &        ' chain',ichain,' chain group',igr,'. Increase MAXNIT.'
+              return1
+            endif
+          endif
+        else
+c          write (iout,*) "No overlap",i-1
+          back=.false.
+c          nit=0 
+c          i=i+1
+        endif
+
+        ENDDO
+
+        if (.not.back) then
+#ifdef DEBUG
+          write (iout,*) "++++++++++Successful generation",igr,ichain,ii
+#endif
+          ii=ii+1
+          nit=0
+        endif
+        back=.false.
+c        write (iout,*) "ii",ii
+        niter=niter+1
+
+      enddo
+      if (niter.ge.maxgen) then
+        write (iout,'(a,2i7,a,i7,a,i7,a,i7)') 
+     & 'Too many trials in conformation generation',niter,maxgen,
+     & ' chain group',igr,' chain',ichain,' residue',ii
+        return1
+      endif
+
+      ENDDO
+
+      do i=2,nres
+        if (itype(i).eq.ntyp1) then
+          do j=1,3
+            dc(j,i-1)=c(j,i)-c(j,i-1)
+          enddo
+        endif
+      enddo
+
+      return
+      end
+c-------------------------------------------------------------------------
+      logical function overlap_mchain(i,ii,ichain,igr)
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      double precision redfac /0.5D0/,redfacp /0.8d0/,redfacscp /0.8d0/
+      integer i,ii,ichain,igr,j,jj,jchain,jgr,k,iti,itj,iteli,itelj
+      double precision rcomp
+      double precision dist
+      logical lprn /.false./
+      overlap_mchain=.false.
+      iti=iabs(itype(i))
+      if (iti.gt.ntyp) return
+C Check for SC-SC overlaps.
+cd    print *,'nnt=',nnt,' nct=',nct
+#ifdef DEBUG
+      write (iout,*) "overlap_mchain i",i," ii",ii," ichain",ichain,
+     &   " igr",igr," iti",iti
+#endif
+      do j=nnt,nct
+        if (itype(j).eq.ntyp1) cycle
+        jchain=ireschain(j)
+c        write (iout,*) "overlap_mchain j",j," jj",jj," jchain",jchain,
+c     &    " itype",itype(j)
+        jgr=mapchain(jchain)
+        jj=j-chain_border1(1,jchain)+1
+        itj=iabs(itype(j))
+c        write (iout,*) "jgr",jgr
+        if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-1
+     &    .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
+        if (j.lt.i-1 .or. ipot.ne.4) then
+          rcomp=sigmaii(iti,itj)
+        else 
+          rcomp=sigma(iti,itj)
+        endif
+cd      print *,'j=',j
+        if (dist(nres+i,nres+j).lt.redfac*rcomp) then
+          overlap_mchain=.true.
+          if (lprn) write(iout,*)'overlap_mchain, SC-SC: i=',i,' j=',j,
+     &     ' ichain',ichain,' jchain',jchain,
+     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
+     &     rcomp*redfac
+        return
+        endif
+      enddo
+#ifdef CHUJ
+C Check for overlaps between the added peptide group and the preceding
+C SCs.
+      iteli=itel(i)
+      if (iteli.gt.0) then
+      do j=1,3
+        c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+      enddo
+      do j=nnt,nct
+        itj=iabs(itype(j))
+        if (itj.eq.ntyp1) cycle
+        jchain=ireschain(j)
+        jgr=mapchain(jchain)
+        jj=j-chain_border1(1,jchain)+1
+        if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2
+     &    .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
+cd      print *,'overlap_mchain, p-Sc: i=',i,' j=',j,
+cd   &         ' dist=',dist(nres+j,maxres2+1)
+        if (dist(j,maxres2+1).lt.4.0D0*redfacscp) then
+          if (lprn) write (iout,*) 'overlap_mchain, p-SC: i=',i,' j=',j,
+     &     ' ichain',ichain,' jchain',jchain,
+     &     ' dist=',dist(nres+j,maxres2+1),' rcomp=',
+     &    4.0d0*redfac 
+          overlap_mchain=.true.
+          return
+        endif
+      enddo
+      endif
+C Check for p-p overlaps
+#endif
+      iteli=itel(i)
+      if (iteli.eq.0) return
+      do j=1,3
+        c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+      enddo
+c      do j=nnt,i-2
+      do j=nnt,nct
+        itelj=itel(j)
+        if (itelj.eq.0) cycle
+        do k=1,3
+          c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
+        enddo
+        jchain=ireschain(j)
+        jgr=mapchain(jchain)
+c        if(iteli.ne.0.and.itelj.ne.0)then
+c        write (iout,*) i,j,dist(maxres2+1,maxres2+2),rpp(iteli,itelj)
+        jj=j-chain_border1(1,jchain)+1
+c        write (iout,*) "jgr",jgr
+        if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2
+     &    .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
+#ifdef DEBUG
+          write (iout,*)'overlap_mchain, p-p: i=',i,' j=',j,' igr',igr,
+     &   ' jgr',jgr,' ichain',ichain,' jchain',jchain,
+     &   ' dist=',dist(maxres2+1,maxres2+2)
+#endif
+        if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfacp) then
+          if (lprn) write (iout,*) 'overlap_mchain, p-p: i=',i,' j=',j,
+     &     ' ichain',ichain,' jchain',jchain,
+     &     ' dist=',dist(maxres2+1,maxres2+2),' rcomp=',
+     &    rpp(iteli,itelj)*redfacp
+          overlap_mchain=.true.
+          return
+        endif
+c        endif
+      enddo
+      return
+      end
+