1 subroutine gen_rand_conf_mchain(nstart0,*)
2 C Generate random conformation or chain cut and regrowth.
8 include 'COMMON.INTERACT'
9 include 'COMMON.IOUNITS'
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
18 double precision gen_theta,gen_phi,dist,ran_number,scalar
19 c write (iout,*) 'gen_rand_conf_mchain: maxgen=',maxgen
22 c write (iout,*) 'Gen_Rand_conf_mchain: nstart=',nstart,
23 c & " nchain_group",nchain_group
27 DO IEQUI=1,NEQUIV(IGR)
29 ichain=iequiv(iequi,igr)
31 i=chain_border1(1,ichain)+nstart-1
34 c(j,i)=ran_number(-15.0d0,15.0d0)
41 dc_norm(j,i)=ran_number(-1.0d0,1.0d0)
43 aux=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i)))
45 dc_norm(j,i)=dc_norm(j,i)/aux
47 if (itype(i).eq.ntyp1) then
49 dc(j,i)=1.9d0*dc_norm(j,i)
53 dc(j,i)=3.8d0*dc_norm(j,i)
57 c(j,i+1)=c(j,i)+dc(j,i)
60 c if (nstart.lt.5) then
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))
69 do while (fail.and.nsi.le.maxsi)
70 call gen_side(it1,theta(i+2),alph(i+1),omeg(i+1),fail)
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.'
80 call orig_frame_chain(i)
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
88 write (iout,*) "dc_norm(:",i-1,") and prod"
90 write (iout,*) j,dc_norm(j,i-1),(prod(j,k,i-1),k=1,3)
105 nnres=chain_border1(2,iequiv(1,igr))-
106 & chain_border1(1,iequiv(1,igr))+1
108 write (iout,*) "chain group",igr," chains",
109 & (iequiv(j,igr),j=1,nequiv(igr))
110 write (iout,*) "ii",ii," nnres",nnres
112 do while (ii.le. nnres .and. niter.lt.maxgen)
115 i=ii-1+chain_border1(1,ichain)
117 write (iout,*) "ii",ii," nnres",nnres," ichain",ichain," i",i,
118 & "niter",niter," back",back," nstart",nstart
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...'
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)
139 write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,
140 & ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen
142 phi(i)=gen_phi(i,it2,it1)
143 phi(i+1)=gen_phi(i+1,it1,it)
145 write (iout,*) "phi",i,phi(i)," phi",i+1,phi(i+1)
147 do iequi=2,nequiv(igr)
148 iii=ii+chain_border1(1,iequiv(iequi,igr))
153 phi(i)=gen_phi(i+1,it2,it1)
155 write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it
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))
162 theta(iii-1)=theta(i-1)
164 if (it2.ne.10 .and. it2.ne.ntyp1) then
167 do while (fail.and.nsi.le.maxsi)
168 call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
171 do iequi=2,nequiv(igr)
172 iii=ii-3+chain_border1(1,iequiv(iequi,igr))
177 write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
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.'
187 c call locate_next_res(i-1)
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
194 theta(i)=ran_number(1.326d0,2.548d0)
197 write (iout,*) "ii",ii," i",i," it1",it1,
198 & " theta",theta(i)," phi",
201 if (it1.ne.10 .and. it1.ne.ntyp1) then
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)
210 write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
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.'
222 do iequi=2,nequiv(igr)
223 iii=ii-1+chain_border1(1,iequiv(iequi,igr))
226 alph(iii-1)=alph(i-1)
227 omeg(iii-1)=omeg(i-1)
230 DO IEQUI=1,NEQUIV(IGR)
232 ichain=iequiv(iequi,igr)
233 i=ii-1+chain_border1(1,ichain)
235 if (back) call locate_next_res(i-1)
237 call locate_next_res(i)
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)
244 if (overlap_mchain(i-1,ii-1,ichain,igr)) then
246 write (iout,*) "***********overlap",i-1," nit",nit
248 if (nit.lt.maxnit) then
254 write (iout,*) "***********overlap maxnit exceeded",nit
260 c write (iout,*) "ii",ii
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.'
270 c write (iout,*) "No overlap",i-1
280 write (iout,*) "++++++++++Successful generation",igr,ichain,ii
286 c write (iout,*) "ii",ii
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
300 if (itype(i).eq.ntyp1) then
302 dc(j,i-1)=c(j,i)-c(j,i-1)
309 c-------------------------------------------------------------------------
310 logical function overlap_mchain(i,ii,ichain,igr)
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.
324 if (iti.gt.ntyp) return
325 C Check for SC-SC overlaps.
326 cd print *,'nnt=',nnt,' nct=',nct
328 write (iout,*) "overlap_mchain i",i," ii",ii," ichain",ichain,
329 & " igr",igr," iti",iti
332 if (itype(j).eq.ntyp1) cycle
334 c write (iout,*) "overlap_mchain j",j," jj",jj," jchain",jchain,
335 c & " itype",itype(j)
337 jj=j-chain_border1(1,jchain)+1
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)
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=',
358 C Check for overlaps between the added peptide group and the preceding
363 c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
367 if (itj.eq.ntyp1) cycle
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=',
380 overlap_mchain=.true.
385 C Check for p-p overlaps
388 if (iteli.eq.0) return
390 c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
395 if (itelj.eq.0) cycle
397 c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
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
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)
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.