1 subroutine gen_rand_conf(nstart,*)
2 C Generate random conformation or chain cut and regrowth.
3 implicit real*8 (a-h,o-z)
8 include 'COMMON.INTERACT'
9 include 'COMMON.IOUNITS'
12 include 'COMMON.CONTROL'
13 logical overlap,back,fail
14 cd print *,' CG Processor',me,' maxgen=',maxgen
15 c write(iout,*) "czy kurwa wogole wchodze"
17 cd write (iout,*) 'Gen_Rand_conf: nstart=',nstart
20 phi(4)=gen_phi(4,iabs(itype(2)),abs(itype(3)))
21 c write(iout,*)'phi(4)=',rad2deg*phi(4)
22 ichir1=isign(1,itype(1))
23 ichir2=isign(1,itype(3))
24 if (nstart.lt.3) theta(3)=gen_theta(itype(2),ichir1,ichir2,
26 write(iout,*)'theta(3)=',rad2deg*theta(3)
30 do while (fail.and.nsi.le.maxsi)
31 call gen_side(it1,theta(3),alph(2),omeg(2),fail)
34 if (nsi.gt.maxsi) return1
49 do while (i.le.nres .and. niter.lt.maxgen)
52 write (iout,'(/80(1h*)/2a/80(1h*))')
53 & 'Generation procedure went down to ',
54 & 'chain beginning. Cannot continue...'
55 write (*,'(/80(1h*)/2a/80(1h*))')
56 & 'Generation procedure went down to ',
57 & 'chain beginning. Cannot continue...'
64 ichir3=isign(1,itype(i))
65 ichir2=isign(1,itype(i-1))
66 ichir0=isign(1,itype(i-3))
67 ichir1=isign(1,itype(i-2))
68 c print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
69 c & ' nit=',nit,' niter=',niter,' maxgen=',maxgen
70 phi(i+1)=gen_phi(i+1,it1,it)
72 phi(i)=gen_phi(i+1,it2,it1)
73 print *,'phi(',i,')=',phi(i)
74 theta(i-1)=gen_theta(it2,ichir0,ichir2,phi(i-1),phi(i))
78 do while (fail.and.nsi.le.maxsi)
79 call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
82 if (nsi.gt.maxsi) return1
84 call locate_next_res(i-1)
86 theta(i)=gen_theta(it1,ichir1,ichir3,phi(i),phi(i+1))
90 do while (fail.and.nsi.le.maxsi)
91 call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
94 if (nsi.gt.maxsi) return1
96 call locate_next_res(i)
97 if (overlap(i-1)) then
98 if (nit.lt.maxnit) then
108 & 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
110 & 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
121 if (niter.ge.maxgen) then
122 write (iout,'(a,2i5)')
123 & 'Too many trials in conformation generation',niter,maxgen
125 & 'Too many trials in conformation generation',niter,maxgen
130 c(j,nres+nres)=c(j,nres)
134 c-------------------------------------------------------------------------
135 logical function overlap(i)
136 implicit real*8 (a-h,o-z)
138 include 'COMMON.CHAIN'
139 include 'COMMON.INTERACT'
140 include 'COMMON.FFIELD'
144 if (iti.gt.ntyp) return
145 C Check for SC-SC overlaps.
146 cd print *,'nnt=',nnt,' nct=',nct
149 if (j.lt.i-1 .or. ipot.ne.4) then
150 rcomp=sigmaii(iti,itj)
155 if (dist(nres+i,nres+j).lt.redfac*rcomp) then
157 c print *,'overlap, SC-SC: i=',i,' j=',j,
158 c & ' dist=',dist(nres+i,nres+j),' rcomp=',
163 C Check for overlaps between the added peptide group and the preceding
167 c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
171 cd print *,'overlap, p-Sc: i=',i,' j=',j,
172 cd & ' dist=',dist(nres+j,maxres2+1)
173 if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
178 C Check for overlaps between the added side chain and the preceding peptide
182 c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1))
184 cd print *,'overlap, SC-p: i=',i,' j=',j,
185 cd & ' dist=',dist(nres+i,maxres2+1)
186 if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then
191 C Check for p-p overlaps
193 c(j,maxres2+2)=0.5D0*(c(j,i)+c(j,i+1))
198 c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
200 cd print *,'overlap, p-p: i=',i,' j=',j,
201 cd & ' dist=',dist(maxres2+1,maxres2+2)
202 if(iteli.ne.0.and.itelj.ne.0)then
203 if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then
211 c--------------------------------------------------------------------------
212 double precision function gen_phi(i,it1,it2)
213 implicit real*8 (a-h,o-z)
216 include 'COMMON.BOUNDS'
217 c gen_phi=ran_number(-pi,pi)
218 C 8/13/98 Generate phi using pre-defined boundaries
219 gen_phi=ran_number(phibound(1,i),phibound(2,i))
222 c---------------------------------------------------------------------------
223 double precision function gen_theta(it,ichir1,ichir2,gama,gama1)
224 implicit real*8 (a-h,o-z)
226 include 'COMMON.LOCAL'
228 double precision y(2),z(2)
229 double precision theta_max,theta_min
230 print *,'gen_theta: it=',it
233 if (dabs(gama).gt.dwapi) then
240 if (dabs(gama1).gt.dwapi) then
247 thet_pred_mean=a0thet(it)
249 thet_pred_mean=thet_pred_mean+athet(k,it,ichir1,ichir2)
250 & *y(k)+bthet(k,it,ichir1,ichir2)*z(k)
254 sig=sig*thet_pred_mean+polthet(j,it)
256 sig=0.5D0/(sig*sig+sigc0(it))
258 &0.5D0*((gthet(2,it)-thet_pred_mean)
260 print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
261 print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
262 theta_temp=binorm(thet_pred_mean,theta0(it),sig
264 if (theta_temp.lt.theta_min) theta_temp=theta_min
265 if (theta_temp.gt.theta_max) theta_temp=theta_max
267 print '(a)','Exiting GENTHETA.'
270 c-------------------------------------------------------------------------
271 subroutine gen_side(it,the,al,om,fail)
272 implicit real*8 (a-h,o-z)
275 include 'COMMON.LOCAL'
276 include 'COMMON.SETUP'
277 include 'COMMON.IOUNITS'
278 double precision MaxBoxLen /10.0D0/
279 double precision Ap_inv(3,3),a(3,3),z(3,maxlob),W1(maxlob),
280 & sumW(0:maxlob),y(2),cm(2),eig(2),box(2,2),work(100),detAp(maxlob)
281 double precision eig_limit /1.0D-8/
282 double precision Big /10.0D0/
283 double precision vec(3,3)
284 logical lprint,fail,lcheck
288 if (the.eq.0.0D0 .or. the.eq.pi) then
290 write (*,'(a,i4,a,i3,a,1pe14.5)')
291 & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
293 cd write (iout,'(a,i3,a,1pe14.5)')
294 cd & 'Error in GenSide: it=',it,' theta=',the
300 nlobit=nlob(iabs(it))
303 print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
304 write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
306 print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
307 write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,
311 zz1=tant-censc(1,i,it)
314 a(k,l)=gaussc(k,l,i,it)
317 detApi=a(2,2)*a(3,3)-a(2,3)**2
318 Ap_inv(2,2)=a(3,3)/detApi
319 Ap_inv(2,3)=-a(2,3)/detApi
320 Ap_inv(3,2)=Ap_inv(2,3)
321 Ap_inv(3,3)=a(2,2)/detApi
323 write (*,'(/a,i2/)') 'Cluster #',i
324 write (*,'(3(1pe14.5),5x,1pe14.5)')
325 & ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
326 write (iout,'(/a,i2/)') 'Cluster #',i
327 write (iout,'(3(1pe14.5),5x,1pe14.5)')
328 & ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
333 W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
337 W1(i)=dexp(bsc(i,iabs(it))-0.5D0*W1i*zz1*zz1)
338 c if (lprint) write(*,'(a,3(1pe15.5)/)')
339 c & 'detAp, W1, anormi',detApi,W1i,anormi
343 zk=zk+zz1*Ap_inv(k,l)*a(l,1)
347 detAp(i)=dsqrt(detApi)
351 print *,'W1:',(w1(i),i=1,nlobit)
352 print *,'detAp:',(detAp(i),i=1,nlobit)
355 print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
357 write (iout,*) 'W1:',(w1(i),i=1,nlobit)
358 write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
361 write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
365 C Writing the distribution just to check the procedure
371 fac=fac+W1(i)/detAp(i)
373 fac=1.0D0/(2.0D0*fac*pi)
374 cd print *,it,'fac=',fac
383 a(j-1,k-1)=gaussc(j,k,i,it)
395 wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
398 wart=wart+W1(i)*dexp(-0.5D0*wykl)
405 c print *,'y',y(1),y(2),' fac=',fac
407 write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
412 c print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
416 C Calculate the CM of the system
423 sumW(i)=sumW(i-1)+W1(i)
428 cm(1)=cm(1)+z(2,j)*W1(j)
429 cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
431 cm(1)=cm(1)/sumW(nlobit)
432 cm(2)=cm(2)/sumW(nlobit)
433 if (cm(1).gt.Big .or. cm(1).lt.-Big .or.
434 & cm(2).gt.Big .or. cm(2).lt.-Big) then
435 cd write (iout,'(a)')
436 cd & 'Unexpected error in GenSide - CM coordinates too large.'
437 cd write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
439 cd & 'Unexpected error in GenSide - CM coordinates too large.'
440 cd write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
444 cd print *,'CM:',cm(1),cm(2)
446 C Find the largest search distance from CM
452 a(j-1,k-1)=gaussc(j,k,i,it)
456 call f02faf('N','U',2,a,3,eig,work,100,ifail)
458 call djacob(2,3,10000,1.0d-10,a,vec,eig)
462 print *,'*************** CG Processor',me
463 print *,'CM:',cm(1),cm(2)
464 write (iout,*) '*************** CG Processor',me
465 write (iout,*) 'CM:',cm(1),cm(2)
466 print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
467 write (iout,'(A,8f10.5)')
468 & 'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
471 if (eig(1).lt.eig_limit) then
473 & 'From Mult_Norm: Eigenvalues of A are too small.'
475 & 'From Mult_Norm: Eigenvalues of A are too small.'
482 radius=radius+pinorm(z(j+1,i)-cm(j))**2
484 radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
485 if (radius.gt.radmax) radmax=radius
487 if (radmax.gt.pi) radmax=pi
489 C Determine the boundaries of the search rectangle.
492 print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
493 print '(a,4(1pe14.4))','radmax: ',radmax
495 box(1,1)=dmax1(cm(1)-radmax,0.0D0)
496 box(2,1)=dmin1(cm(1)+radmax,pi)
497 box(1,2)=cm(2)-radmax
498 box(2,2)=cm(2)+radmax
501 print *,'CG Processor',me,' Array BOX:'
505 print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
506 print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
508 write (iout,*)'CG Processor',me,' Array BOX:'
510 write (iout,*)'Array BOX:'
512 write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
513 write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
515 if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
517 write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
518 write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
520 c write (iout,'(a)') 'Bad sampling box.'
525 which_lobe=ran_number(0.0D0,sumW(nlobit))
526 c print '(a,1pe14.4)','which_lobe=',which_lobe
528 if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
531 c print *,'ilob=',ilob,' nlob=',nlob(it)
535 a(i-1,j-1)=gaussc(i,j,ilob,it)
538 cd print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
539 call mult_norm1(3,2,a,cm,box,y,fail)
543 cd print *,'al=',al,' om=',om
547 c---------------------------------------------------------------------------
548 double precision function ran_number(x1,x2)
549 C Calculate a random real number from the range (x1,x2).
550 implicit real*8 (a-h,o-z)
552 double precision x1,x2,fctor
553 data fctor /2147483647.0D0/
556 include 'COMMON.SETUP'
557 ran_number=x1+(x2-x1)*prng_next(me)
560 ran_number=x1+(x2-x1)*ix/fctor
564 c--------------------------------------------------------------------------
565 integer function iran_num(n1,n2)
566 C Calculate a random integer number from the range (n1,n2).
567 implicit real*8 (a-h,o-z)
570 real fctor /2147483647.0/
573 include 'COMMON.SETUP'
574 ix=n1+(n2-n1+1)*prng_next(me)
580 ix=n1+(n2-n1+1)*(ix/fctor)
586 c--------------------------------------------------------------------------
587 double precision function binorm(x1,x2,sigma1,sigma2,ak)
588 implicit real*8 (a-h,o-z)
589 print *,'Enter BINORM.',x1,x2,sigma1,sigma2,ak
590 alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
591 aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
592 seg=sigma1/(sigma1+ak*sigma2)
593 alen=ran_number(0.0D0,1.0D0)
594 if (alen.lt.seg) then
595 c print *,'przed anorm',x1,sigma1,alowb,aupb
596 c print *, 'anorm',anorm_distr(x1,sigma1,alowb,aupb)
597 binorm=anorm_distr(x1,sigma1,alowb,aupb)
600 c print *,'przed anorm',x2,sigma2,alowb,aupb
601 c print *, 'anorm',anorm_distr(x2,sigma2,alowb,aupb)
602 binorm=anorm_distr(x2,sigma2,alowb,aupb)
604 print '(a)','Exiting BINORM.'
607 c-----------------------------------------------------------------------
608 c double precision function anorm_distr(x,sigma,alowb,aupb)
609 c implicit real*8 (a-h,o-z)
610 c print '(a)','Enter ANORM_DISTR.'
611 c 10 y=ran_number(alowb,aupb)
612 c expon=dexp(-0.5D0*((y-x)/sigma)**2)
613 c ran=ran_number(0.0D0,1.0D0)
614 c if (expon.lt.ran) goto 10
616 c print '(a)','Exiting ANORM_DISTR.'
619 c-----------------------------------------------------------------------
620 double precision function anorm_distr(x,sigma,alowb,aupb)
621 implicit real*8 (a-h,o-z)
622 c to make a normally distributed deviate with zero mean and unit variance
625 real fac,gset,rsq,v1,v2,ran1
629 1 v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
630 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
631 c print *,"anorm: iset",iset," v1",v1," v2",v2," rsq",rsq
633 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
634 fac=sqrt(-2.0d0*log(rsq)/rsq)
642 anorm_distr=x+gaussdev*sigma
645 c------------------------------------------------------------------------
646 subroutine mult_norm(lda,n,a,x,fail)
648 C Generate the vector X whose elements obey the multiple-normal distribution
649 C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
650 C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
651 C eigenvalue of the matrix A is close to 0.
653 implicit double precision (a-h,o-z)
654 double precision a(lda,n),x(n),eig(100),vec(3,3),work(100)
655 double precision eig_limit /1.0D-8/
658 c print '(a)','Enter MULT_NORM.'
660 C Find the smallest eigenvalue of the matrix A.
663 c print '(8f10.5)',(a(i,j),j=1,n)
666 call f02faf('V','U',2,a,lda,eig,work,100,ifail)
668 call djacob(2,lda,10000,1.0d-10,a,vec,eig)
670 c print '(8f10.5)',(eig(i),i=1,n)
673 c print '(8f10.5)',(a(i,j),j=1,n)
675 if (eig(1).lt.eig_limit) then
676 print *,'From Mult_Norm: Eigenvalues of A are too small.'
681 C Generate points following the normal distributions along the principal
682 C axes of the moment matrix. Store in WORK.
685 sigma=1.0D0/dsqrt(eig(i))
687 work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
690 C Transform the vector of normal variables back to the original basis.
701 c------------------------------------------------------------------------
702 subroutine mult_norm1(lda,n,a,z,box,x,fail)
704 C Generate the vector X whose elements obey the multi-gaussian multi-dimensional
705 C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
706 C leading dimension of the moment matrix A, n is the dimension of the
707 C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
708 C smallest eigenvalue of the matrix A is close to 0.
710 implicit real*8 (a-h,o-z)
715 double precision a(lda,n),z(n),x(n),box(n,n)
716 double precision etmp
717 include 'COMMON.IOUNITS'
719 include 'COMMON.SETUP'
723 C Generate points following the normal distributions along the principal
724 C axes of the moment matrix. Store in WORK.
726 cd print *,'CG Processor',me,' entered MultNorm1.'
727 cd print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
729 cd print *,i,box(1,i),box(2,i)
733 if (istep.gt.10000) then
734 c write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
736 c write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
738 c write (iout,*) 'box',box
739 c write (iout,*) 'a',a
740 c write (iout,*) 'z',z
745 x(i)=ran_number(box(1,i),box(2,i))
750 ww=ww+0.5D0*a(i,i)*xi*xi
752 ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
755 dec=ran_number(0.0D0,1.0D0)
756 c print *,(x(i),i=1,n),ww,dexp(-ww),dec
757 crc if (dec.gt.dexp(-ww)) goto 10
763 if (dec.gt.etmp) goto 10
764 cd print *,'CG Processor',me,' exitting MultNorm1.'
768 crc--------------------------------------
769 subroutine overlap_sc(scfail)
770 c Internal and cartesian coordinates must be consistent as input,
771 c and will be up-to-date on return.
772 c At the end of this procedure, scfail is true if there are
773 c overlapping residues left, or false otherwise (success)
774 implicit real*8 (a-h,o-z)
776 include 'COMMON.CHAIN'
777 include 'COMMON.INTERACT'
778 include 'COMMON.FFIELD'
780 include 'COMMON.SBRIDGE'
781 include 'COMMON.IOUNITS'
782 logical had_overlaps,fail,scfail
783 integer ioverlap(maxres),ioverlap_last
786 call overlap_sc_list(ioverlap,ioverlap_last)
787 if (ioverlap_last.gt.0) then
788 write (iout,*) '#OVERLAPing residues ',ioverlap_last
789 write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
795 if (ioverlap_last.eq.0) exit
797 do ires=1,ioverlap_last
803 do while (fail.and.nsi.le.maxsi)
804 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
812 call overlap_sc_list(ioverlap,ioverlap_last)
813 c write (iout,*) 'Overlaping residues ',ioverlap_last,
814 c & (ioverlap(j),j=1,ioverlap_last)
817 if (k.le.1000.and.ioverlap_last.eq.0) then
819 if (had_overlaps) then
820 write (iout,*) '#OVERLAPing all corrected after ',k,
821 & ' random generation'
825 write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
826 write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
832 write (iout,'(a30,i5,a12,i4)')
833 & '#OVERLAP FAIL in gen_side after',maxsi,
839 subroutine overlap_sc_list(ioverlap,ioverlap_last)
840 implicit real*8 (a-h,o-z)
843 include 'COMMON.LOCAL'
844 include 'COMMON.IOUNITS'
845 include 'COMMON.CHAIN'
846 include 'COMMON.INTERACT'
847 include 'COMMON.FFIELD'
849 include 'COMMON.CALC'
851 integer ioverlap(maxres),ioverlap_last
855 C Check for SC-SC overlaps and mark residues
856 c print *,'>>overlap_sc nnt=',nnt,' nct=',nct
860 itypi1=abs(itype(i+1))
864 dxi=dc_norm(1,nres+i)
865 dyi=dc_norm(2,nres+i)
866 dzi=dc_norm(3,nres+i)
867 dsci_inv=dsc_inv(itypi)
870 do j=istart(i,iint),iend(i,iint)
873 dscj_inv=dsc_inv(itypj)
874 sig0ij=sigma(itypi,itypj)
875 chi1=chi(itypi,itypj)
876 chi2=chi(itypj,itypi)
883 alf12=0.5D0*(alf1+alf2)
885 rcomp=sigmaii(itypi,itypj)
887 rcomp=sigma(itypi,itypj)
889 c print '(2(a3,2i3),a3,2f10.5)',
890 c & ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
895 dxj=dc_norm(1,nres+j)
896 dyj=dc_norm(2,nres+j)
897 dzj=dc_norm(3,nres+j)
898 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
902 sig=sig0ij*dsqrt(sigsq)
903 rij_shift=1.0D0/rij-sig+sig0ij
905 ct if ( 1.0/rij .lt. redfac*rcomp .or.
906 ct & rij_shift.le.0.0D0 ) then
907 if ( rij_shift.le.0.0D0 ) then
908 cd write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
909 cd & 'overlap SC-SC: i=',i,' j=',j,
910 cd & ' dist=',dist(nres+i,nres+j),' rcomp=',
911 cd & rcomp,1.0/rij,rij_shift
912 ioverlap_last=ioverlap_last+1
913 ioverlap(ioverlap_last)=i
914 do k=1,ioverlap_last-1
915 if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
917 ioverlap_last=ioverlap_last+1
918 ioverlap(ioverlap_last)=j
919 do k=1,ioverlap_last-1
920 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1