--- /dev/null
+ 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
+
projects / unres.git / blobdiff