source/unres/src-HCD-5D/gen_rand_conf_mchain.F

   1       subroutine gen_rand_conf_mchain(nstart0,*)
   2 C Generate random conformation or chain cut and regrowth.
   3       implicit none
   4       include 'DIMENSIONS'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.LOCAL'
   7       include 'COMMON.VAR'
   8       include 'COMMON.INTERACT'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.MCM'
  11       include 'COMMON.GEO'
  12       include 'COMMON.CONTROL'
  13       logical overlap_mchain,back,fail
  14       integer nstart,nstart0
  15       integer i,ii,iii,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit
  16       integer igr,iequi,ichain,nnres
  17       double precision aux
  18       double precision gen_theta,gen_phi,dist,ran_number,scalar
  19 c      write (iout,*)  'gen_rand_conf_mchain: maxgen=',maxgen
  20       nstart=nstart0
  21       maxsi=100
  22 c      write (iout,*) 'Gen_Rand_conf_mchain: nstart=',nstart,
  23 c     &   " nchain_group",nchain_group
  24
  25       DO IGR=1,NCHAIN_GROUP
  26
  27       DO IEQUI=1,NEQUIV(IGR)
  28
  29       ichain=iequiv(iequi,igr)
  30
  31       i=chain_border1(1,ichain)+nstart-1
  32       if (nstart.eq.1) then
  33         do j=1,3
  34           c(j,i)=ran_number(-15.0d0,15.0d0)
  35           dc(j,i-1)=c(j,i)
  36         enddo
  37       endif
  38       if (nstart.le.2) then
  39
  40       do j=1,3
  41         dc_norm(j,i)=ran_number(-1.0d0,1.0d0)
  42       enddo
  43       aux=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i)))
  44       do j=1,3
  45         dc_norm(j,i)=dc_norm(j,i)/aux
  46       enddo
  47       if (itype(i).eq.ntyp1) then
  48         do j=1,3
  49           dc(j,i)=1.9d0*dc_norm(j,i)
  50         enddo
  51       else
  52         do j=1,3
  53           dc(j,i)=3.8d0*dc_norm(j,i)
  54         enddo
  55       endif
  56       do j=1,3
  57         c(j,i+1)=c(j,i)+dc(j,i)
  58       enddo
  59       endif
  60 c      if (nstart.lt.5) then
  61       if (nstart.le.2) then
  62       it1=iabs(itype(i+2))
  63       phi(i+3)=gen_phi(i+3,iabs(itype(i+1)),iabs(itype(i+2)))
  64 c       write(iout,*)'phi(4)=',rad2deg*phi(4)
  65       theta(i+2)=gen_theta(iabs(itype(i+2)),pi,phi(i+3))
  66       if (it1.ne.10) then
  67         nsi=0
  68         fail=.true.
  69         do while (fail.and.nsi.le.maxsi)
  70           call gen_side(it1,theta(i+2),alph(i+1),omeg(i+1),fail)
  71           nsi=nsi+1
  72         enddo
  73         if (nsi.gt.maxsi) then
  74           write (iout,'(a,i7,a,i7,a,i7,a,i7)')
  75      &    'Problem in SC rotamer generation, residue',ii,
  76      &    ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
  77           return1
  78         endif
  79       endif ! it1.ne.10
  80       call orig_frame_chain(i)
  81       else
  82 c      call orig_frame_chain(i-1)
  83 c      write (iout,*) "calling refsys",i
  84       call refsys(i,i-1,i-2,prod(1,1,i-1),
  85      &  prod(1,2,i-1),prod(1,3,i-1),fail)
  86 c      write (iout,*) "after refsys",i
  87 #ifdef DEBUG
  88       write (iout,*) "dc_norm(:",i-1,") and prod"
  89       do j=1,3
  90         write (iout,*) j,dc_norm(j,i-1),(prod(j,k,i-1),k=1,3)
  91       enddo
  92 #endif
  93       endif
  94
  95       ENDDO
  96
  97       nstart=nstart+1
  98       ii=nstart
  99
 100       maxnit=5000
 101
 102       nit=0
 103       niter=0
 104       back=.false.
 105       nnres=chain_border1(2,iequiv(1,igr))-
 106      &  chain_border1(1,iequiv(1,igr))+1
 107 #ifdef DEBUG
 108       write (iout,*) "chain group",igr," chains",
 109      &    (iequiv(j,igr),j=1,nequiv(igr))
 110       write (iout,*) "ii",ii," nnres",nnres
 111 #endif
 112       do while (ii.le. nnres .and. niter.lt.maxgen)
 113
 114         ichain=iequiv(1,igr)
 115         i=ii-1+chain_border1(1,ichain)
 116 #ifdef DEBUG
 117         write (iout,*) "ii",ii," nnres",nnres," ichain",ichain," i",i,
 118      &    "niter",niter," back",back," nstart",nstart
 119 #endif
 120         if (ii.lt.nstart) then
 121 c          if(iprint.gt.1) then
 122           write (iout,'(/80(1h*)/2a/80(1h*))')
 123      &          'Generation procedure went down to ',
 124      &          'chain beginning. Cannot continue...'
 125           write (*,'(/80(1h*)/2a/80(1h*))')
 126      &          'Generation procedure went down to ',
 127      &          'chain beginning. Cannot continue...'
 128 c          endif
 129           return1
 130         endif
 131         it1=iabs(itype(i-1))
 132         it2=iabs(itype(i-2))
 133         it=iabs(itype(i))
 134         if (it.eq.ntyp1 .and. it1.eq.ntyp1) then
 135           vbld(i)=ran_number(3.8d0,10.0d0)
 136           vbld_inv(i)=1.0d0/vbld(i)
 137         endif
 138 #ifdef DEBUG
 139         write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,
 140      &    ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen
 141 #endif
 142         phi(i)=gen_phi(i,it2,it1)
 143         phi(i+1)=gen_phi(i+1,it1,it)
 144 #ifdef DEBUG
 145         write (iout,*) "phi",i,phi(i)," phi",i+1,phi(i+1)
 146 #endif
 147         do iequi=2,nequiv(igr)
 148           iii=ii+chain_border1(1,iequiv(iequi,igr))
 149           phi(iii)=phi(i+1)
 150         enddo
 151 #ifdef CHUJ
 152         if (back) then
 153           phi(i)=gen_phi(i+1,it2,it1)
 154 #ifdef DEBUG
 155           write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it
 156 #endif
 157           theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
 158           if (theta(i-1).gt.2.68780478d0) theta(i-1)=2.68780478d0
 159           do iequi=2,nequiv(igr)
 160             iii=ii-1+chain_border1(1,iequiv(iequi,igr))
 161             phi(iii)=phi(i+1)
 162             theta(iii-1)=theta(i-1)
 163           enddo
 164           if (it2.ne.10 .and. it2.ne.ntyp1) then
 165             nsi=0
 166             fail=.true.
 167             do while (fail.and.nsi.le.maxsi)
 168               call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
 169               nsi=nsi+1
 170             enddo
 171             do iequi=2,nequiv(igr)
 172               iii=ii-3+chain_border1(1,iequiv(iequi,igr))
 173               alph(iii)=alph(i-2)
 174               omeg(iii)=omeg(i-2)
 175             enddo
 176 #ifdef DEBUG
 177             write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
 178 #endif
 179             if (nsi.gt.maxsi) then
 180               write (iout,'(a,i7,a,i7,a,i7,a,i7)')
 181      &        'Problem in SC rotamer generation, residue',ii,
 182      &        ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
 183               return1
 184             endif
 185
 186           endif
 187 c          call locate_next_res(i-1)
 188         endif
 189 #endif
 190         if (it1.ne.ntyp1) then
 191           theta(i)=gen_theta(it1,phi(i),phi(i+1))
 192           if (theta(i).gt.2.68780478d0) theta(i)=2.68780478d0
 193         else
 194           theta(i)=ran_number(1.326d0,2.548d0)
 195         endif
 196 #ifdef DEBUG
 197         write (iout,*) "ii",ii," i",i," it1",it1,
 198      &   " theta",theta(i)," phi",
 199      &    phi(i),phi(i+1)
 200 #endif
 201         if (it1.ne.10 .and. it1.ne.ntyp1) then
 202         nsi=0
 203         fail=.true.
 204         do while (fail.and.nsi.le.maxsi)
 205 c          theta(i)=gen_theta(it1,phi(i),phi(i+1))
 206           call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
 207           nsi=nsi+1
 208         enddo
 209 #ifdef DEBUG
 210         write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
 211 #endif
 212 c        write (iout,*) "After forward SC generation:",nsi,maxsi
 213         if (nsi.gt.maxsi) then
 214           write (iout,'(a,i7,a,i7,a,i7,a,i7)')
 215      &    'Problem in SC rotamer generation, residue',ii,
 216      &    ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
 217           return1
 218         endif
 219
 220         endif
 221
 222         do iequi=2,nequiv(igr)
 223           iii=ii-1+chain_border1(1,iequiv(iequi,igr))
 224           theta(iii)=theta(i)
 225           phi(iii)=phi(i)
 226           alph(iii-1)=alph(i-1)
 227           omeg(iii-1)=omeg(i-1)
 228         enddo
 229
 230         DO IEQUI=1,NEQUIV(IGR)
 231
 232         ichain=iequiv(iequi,igr)
 233         i=ii-1+chain_border1(1,ichain)
 234 #ifdef CHUJ
 235         if (back) call locate_next_res(i-1)
 236 #endif
 237         call locate_next_res(i)
 238 #ifdef DEBUG
 239         write (iout,*) "i",i," ii",ii," ichain",ichain
 240         write (iout,*) theta(i)*rad2deg,phi(i)*rad2deg,
 241      &    alph(i-1)*rad2deg,omeg(i-1)*rad2deg
 242         write (iout,*) (c(j,i),j=1,3),(c(j,i+nres-1),j=1,3)
 243 #endif
 244         if (overlap_mchain(i-1,ii-1,ichain,igr)) then
 245 #ifdef DEBUG
 246           write (iout,*) "***********overlap",i-1," nit",nit
 247 #endif
 248           if (nit.lt.maxnit) then
 249             back=.true.
 250             nit=nit+1
 251             exit
 252           else
 253 #ifdef DEBUG
 254             write (iout,*) "***********overlap maxnit exceeded",nit
 255 #endif
 256             nit=0
 257             if (ii.gt.3) then
 258               back=.true.
 259               ii=ii-1
 260 c              write (iout,*) "ii",ii
 261               exit
 262             else
 263               write (iout,'(a,i7,a,i7,a,i7,a,i7)')
 264      &        'Cannot generate non-overlaping conformation, residue',ii,
 265      &        ' chain',ichain,' chain group',igr,'. Increase MAXNIT.'
 266               return1
 267             endif
 268           endif
 269         else
 270 c          write (iout,*) "No overlap",i-1
 271           back=.false.
 272 c          nit=0
 273 c          i=i+1
 274         endif
 275
 276         ENDDO
 277
 278         if (.not.back) then
 279 #ifdef DEBUG
 280           write (iout,*) "++++++++++Successful generation",igr,ichain,ii
 281 #endif
 282           ii=ii+1
 283           nit=0
 284         endif
 285         back=.false.
 286 c        write (iout,*) "ii",ii
 287         niter=niter+1
 288
 289       enddo
 290       if (niter.ge.maxgen) then
 291         write (iout,'(a,2i7,a,i7,a,i7,a,i7)')
 292      & 'Too many trials in conformation generation',niter,maxgen,
 293      & ' chain group',igr,' chain',ichain,' residue',ii
 294         return1
 295       endif
 296
 297       ENDDO
 298
 299       do i=2,nres
 300         if (itype(i).eq.ntyp1) then
 301           do j=1,3
 302             dc(j,i-1)=c(j,i)-c(j,i-1)
 303           enddo
 304         endif
 305       enddo
 306
 307       return
 308       end
 309 c-------------------------------------------------------------------------
 310       logical function overlap_mchain(i,ii,ichain,igr)
 311       implicit none
 312       include 'DIMENSIONS'
 313       include 'COMMON.IOUNITS'
 314       include 'COMMON.CHAIN'
 315       include 'COMMON.INTERACT'
 316       include 'COMMON.FFIELD'
 317       double precision redfac /0.5D0/,redfacp /0.8d0/,redfacscp /0.8d0/
 318       integer i,ii,ichain,igr,j,jj,jchain,jgr,k,iti,itj,iteli,itelj
 319       double precision rcomp
 320       double precision dist
 321       logical lprn /.false./
 322       overlap_mchain=.false.
 323       iti=iabs(itype(i))
 324       if (iti.gt.ntyp) return
 325 C Check for SC-SC overlaps.
 326 cd    print *,'nnt=',nnt,' nct=',nct
 327 #ifdef DEBUG
 328       write (iout,*) "overlap_mchain i",i," ii",ii," ichain",ichain,
 329      &   " igr",igr," iti",iti
 330 #endif
 331       do j=nnt,nct
 332         if (itype(j).eq.ntyp1) cycle
 333         jchain=ireschain(j)
 334 c        write (iout,*) "overlap_mchain j",j," jj",jj," jchain",jchain,
 335 c     &    " itype",itype(j)
 336         jgr=mapchain(jchain)
 337         jj=j-chain_border1(1,jchain)+1
 338         itj=iabs(itype(j))
 339 c        write (iout,*) "jgr",jgr
 340         if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-1
 341      &    .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
 342         if (j.lt.i-1 .or. ipot.ne.4) then
 343           rcomp=sigmaii(iti,itj)
 344         else
 345           rcomp=sigma(iti,itj)
 346         endif
 347 cd      print *,'j=',j
 348         if (dist(nres+i,nres+j).lt.redfac*rcomp) then
 349           overlap_mchain=.true.
 350           if (lprn) write(iout,*)'overlap_mchain, SC-SC: i=',i,' j=',j,
 351      &     ' ichain',ichain,' jchain',jchain,
 352      &     ' dist=',dist(nres+i,nres+j),' rcomp=',
 353      &     rcomp*redfac
 354         return
 355         endif
 356       enddo
 357 #ifdef CHUJ
 358 C Check for overlaps between the added peptide group and the preceding
 359 C SCs.
 360       iteli=itel(i)
 361       if (iteli.gt.0) then
 362       do j=1,3
 363         c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
 364       enddo
 365       do j=nnt,nct
 366         itj=iabs(itype(j))
 367         if (itj.eq.ntyp1) cycle
 368         jchain=ireschain(j)
 369         jgr=mapchain(jchain)
 370         jj=j-chain_border1(1,jchain)+1
 371         if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2
 372      &    .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
 373 cd      print *,'overlap_mchain, p-Sc: i=',i,' j=',j,
 374 cd   &         ' dist=',dist(nres+j,maxres2+1)
 375         if (dist(j,maxres2+1).lt.4.0D0*redfacscp) then
 376           if (lprn) write (iout,*) 'overlap_mchain, p-SC: i=',i,' j=',j,
 377      &     ' ichain',ichain,' jchain',jchain,
 378      &     ' dist=',dist(nres+j,maxres2+1),' rcomp=',
 379      &    4.0d0*redfac
 380           overlap_mchain=.true.
 381           return
 382         endif
 383       enddo
 384       endif
 385 C Check for p-p overlaps
 386 #endif
 387       iteli=itel(i)
 388       if (iteli.eq.0) return
 389       do j=1,3
 390         c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
 391       enddo
 392 c      do j=nnt,i-2
 393       do j=nnt,nct
 394         itelj=itel(j)
 395         if (itelj.eq.0) cycle
 396         do k=1,3
 397           c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
 398         enddo
 399         jchain=ireschain(j)
 400         jgr=mapchain(jchain)
 401 c        if(iteli.ne.0.and.itelj.ne.0)then
 402 c        write (iout,*) i,j,dist(maxres2+1,maxres2+2),rpp(iteli,itelj)
 403         jj=j-chain_border1(1,jchain)+1
 404 c        write (iout,*) "jgr",jgr
 405         if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2
 406      &    .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
 407 #ifdef DEBUG
 408           write (iout,*)'overlap_mchain, p-p: i=',i,' j=',j,' igr',igr,
 409      &   ' jgr',jgr,' ichain',ichain,' jchain',jchain,
 410      &   ' dist=',dist(maxres2+1,maxres2+2)
 411 #endif
 412         if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfacp) then
 413           if (lprn) write (iout,*) 'overlap_mchain, p-p: i=',i,' j=',j,
 414      &     ' ichain',ichain,' jchain',jchain,
 415      &     ' dist=',dist(maxres2+1,maxres2+2),' rcomp=',
 416      &    rpp(iteli,itelj)*redfacp
 417           overlap_mchain=.true.
 418           return
 419         endif
 420 c        endif
 421       enddo
 422       return
 423       end
 424