subroutine gen_rand_conf(nstart,*) C Generate random conformation or chain cut and regrowth. implicit real*8 (a-h,o-z) 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,back,fail cd print *,' CG Processor',me,' maxgen=',maxgen maxsi=100 cd write (iout,*) 'Gen_Rand_conf: nstart=',nstart if (nstart.lt.5) then it1=iabs(itype(2)) phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3))) c write(iout,*)'phi(4)=',rad2deg*phi(4) if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4)) c write(iout,*)'theta(3)=',rad2deg*theta(3) if (it1.ne.10) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) call gen_side(it1,theta(3),alph(2),omeg(2),fail) nsi=nsi+1 enddo if (nsi.gt.maxsi) return1 endif ! it1.ne.10 call orig_frame i=4 nstart=4 else i=nstart nstart=max0(i,4) endif maxnit=0 nit=0 niter=0 back=.false. do while (i.le.nres .and. niter.lt.maxgen) if (i.lt.nstart) then 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...' endif return1 endif it1=iabs(itype(i-1)) it2=iabs(itype(i-2)) it=iabs(itype(i)) c print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2, c & ' nit=',nit,' niter=',niter,' maxgen=',maxgen phi(i+1)=gen_phi(i+1,it1,it) if (back) then phi(i)=gen_phi(i+1,it2,it1) c print *,'phi(',i,')=',phi(i) theta(i-1)=gen_theta(it2,phi(i-1),phi(i)) if (it2.ne.10) 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 if (nsi.gt.maxsi) return1 endif call locate_next_res(i-1) endif theta(i)=gen_theta(it1,phi(i),phi(i+1)) if (it1.ne.10) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail) nsi=nsi+1 enddo if (nsi.gt.maxsi) return1 endif call locate_next_res(i) if (overlap(i-1)) then if (nit.lt.maxnit) then back=.true. nit=nit+1 else nit=0 if (i.gt.3) then back=.true. i=i-1 else write (iout,'(a)') & 'Cannot generate non-overlaping conformation. Increase MAXNIT.' write (*,'(a)') & 'Cannot generate non-overlaping conformation. Increase MAXNIT.' return1 endif endif else back=.false. nit=0 i=i+1 endif niter=niter+1 enddo if (niter.ge.maxgen) then write (iout,'(a,2i5)') & 'Too many trials in conformation generation',niter,maxgen write (*,'(a,2i5)') & 'Too many trials in conformation generation',niter,maxgen return1 endif do j=1,3 c(j,nres+1)=c(j,1) c(j,nres+nres)=c(j,nres) enddo return end c------------------------------------------------------------------------- logical function overlap(i) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.FFIELD' data redfac /0.5D0/ overlap=.false. iti=iabs(itype(i)) if (iti.gt.ntyp) return C Check for SC-SC overlaps. cd print *,'nnt=',nnt,' nct=',nct do j=nnt,i-1 itj=iabs(itype(j)) 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=.true. c print *,'overlap, SC-SC: i=',i,' j=',j, c & ' dist=',dist(nres+i,nres+j),' rcomp=', c & rcomp return endif enddo C Check for overlaps between the added peptide group and the preceding C SCs. iteli=itel(i) do j=1,3 c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1)) enddo do j=nnt,i-2 itj=iabs(itype(j)) cd print *,'overlap, p-Sc: i=',i,' j=',j, cd & ' dist=',dist(nres+j,maxres2+1) if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then overlap=.true. return endif enddo C Check for overlaps between the added side chain and the preceding peptide C groups. do j=1,nnt-2 do k=1,3 c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1)) enddo cd print *,'overlap, SC-p: i=',i,' j=',j, cd & ' dist=',dist(nres+i,maxres2+1) if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then overlap=.true. return endif enddo C Check for p-p overlaps do j=1,3 c(j,maxres2+2)=0.5D0*(c(j,i)+c(j,i+1)) enddo do j=nnt,i-2 itelj=itel(j) do k=1,3 c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1)) enddo cd print *,'overlap, p-p: i=',i,' j=',j, cd & ' dist=',dist(maxres2+1,maxres2+2) if(iteli.ne.0.and.itelj.ne.0)then if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then overlap=.true. return endif endif enddo return end c-------------------------------------------------------------------------- double precision function gen_phi(i,it1,it2) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.BOUNDS' c gen_phi=ran_number(-pi,pi) C 8/13/98 Generate phi using pre-defined boundaries gen_phi=ran_number(phibound(1,i),phibound(2,i)) return end c--------------------------------------------------------------------------- double precision function gen_theta(it,gama,gama1) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.LOCAL' include 'COMMON.GEO' double precision y(2),z(2) double precision theta_max,theta_min c print *,'gen_theta: it=',it theta_min=0.05D0*pi theta_max=0.95D0*pi if (dabs(gama).gt.dwapi) then y(1)=dcos(gama) y(2)=dsin(gama) else y(1)=0.0D0 y(2)=0.0D0 endif if (dabs(gama1).gt.dwapi) then z(1)=dcos(gama1) z(2)=dsin(gama1) else z(1)=0.0D0 z(2)=0.0D0 endif thet_pred_mean=a0thet(it) do k=1,2 thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k) & +bthet(k,it,1,1)*z(k) enddo sig=polthet(3,it) do j=2,0,-1 sig=sig*thet_pred_mean+polthet(j,it) enddo sig=0.5D0/(sig*sig+sigc0(it)) ak=dexp(gthet(1,it)- &0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2) c print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3) c print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak) if (theta_temp.lt.theta_min) theta_temp=theta_min if (theta_temp.gt.theta_max) theta_temp=theta_max gen_theta=theta_temp c print '(a)','Exiting GENTHETA.' return end c------------------------------------------------------------------------- subroutine gen_side(it,the,al,om,fail) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.SETUP' include 'COMMON.IOUNITS' double precision MaxBoxLen /10.0D0/ double precision Ap_inv(3,3),a(3,3),z(3,maxlob),W1(maxlob), & sumW(0:maxlob),y(2),cm(2),eig(2),box(2,2),work(100),detAp(maxlob) double precision eig_limit /1.0D-8/ double precision Big /10.0D0/ double precision vec(3,3) logical lprint,fail,lcheck lcheck=.false. lprint=.false. fail=.false. if (the.eq.0.0D0 .or. the.eq.pi) then #ifdef MPI write (*,'(a,i4,a,i3,a,1pe14.5)') & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the #else cd write (iout,'(a,i3,a,1pe14.5)') cd & 'Error in GenSide: it=',it,' theta=',the #endif fail=.true. return endif tant=dtan(the-pipol) nlobit=nlob(it) if (lprint) then #ifdef MPI print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.' write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.' #endif print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the, & ' tant=',tant endif do i=1,nlobit zz1=tant-censc(1,i,it) do k=1,3 do l=1,3 a(k,l)=gaussc(k,l,i,it) enddo enddo detApi=a(2,2)*a(3,3)-a(2,3)**2 Ap_inv(2,2)=a(3,3)/detApi Ap_inv(2,3)=-a(2,3)/detApi Ap_inv(3,2)=Ap_inv(2,3) Ap_inv(3,3)=a(2,2)/detApi if (lprint) then write (*,'(/a,i2/)') 'Cluster #',i write (*,'(3(1pe14.5),5x,1pe14.5)') & ((a(l,k),l=1,3),censc(k,i,it),k=1,3) write (iout,'(/a,i2/)') 'Cluster #',i write (iout,'(3(1pe14.5),5x,1pe14.5)') & ((a(l,k),l=1,3),censc(k,i,it),k=1,3) endif W1i=0.0D0 do k=2,3 do l=2,3 W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l) enddo enddo W1i=a(1,1)-W1i W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1) c if (lprint) write(*,'(a,3(1pe15.5)/)') c & 'detAp, W1, anormi',detApi,W1i,anormi do k=2,3 zk=censc(k,i,it) do l=2,3 zk=zk+zz1*Ap_inv(k,l)*a(l,1) enddo z(k,i)=zk enddo detAp(i)=dsqrt(detApi) enddo if (lprint) then print *,'W1:',(w1(i),i=1,nlobit) print *,'detAp:',(detAp(i),i=1,nlobit) print *,'Z' do i=1,nlobit print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3) enddo write (iout,*) 'W1:',(w1(i),i=1,nlobit) write (iout,*) 'detAp:',(detAp(i),i=1,nlobit) write (iout,*) 'Z' do i=1,nlobit write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3) enddo endif if (lcheck) then C Writing the distribution just to check the procedure fac=0.0D0 dV=deg2rad**2*10.0D0 sum=0.0D0 sum1=0.0D0 do i=1,nlobit fac=fac+W1(i)/detAp(i) enddo fac=1.0D0/(2.0D0*fac*pi) cd print *,it,'fac=',fac do ial=90,180,2 y(1)=deg2rad*ial do iom=-180,180,5 y(2)=deg2rad*iom wart=0.0D0 do i=1,nlobit do j=2,3 do k=2,3 a(j-1,k-1)=gaussc(j,k,i,it) enddo enddo y2=y(2) do iii=-1,1 y(2)=y2+iii*dwapi wykl=0.0D0 do j=1,2 do k=1,2 wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i)) enddo enddo wart=wart+W1(i)*dexp(-0.5D0*wykl) enddo y(2)=y2 enddo c print *,'y',y(1),y(2),' fac=',fac wart=fac*wart write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart sum=sum+wart sum1=sum1+1.0D0 enddo enddo c print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV return endif C Calculate the CM of the system C do i=1,nlobit W1(i)=W1(i)/detAp(i) enddo sumW(0)=0.0D0 do i=1,nlobit sumW(i)=sumW(i-1)+W1(i) enddo cm(1)=z(2,1)*W1(1) cm(2)=z(3,1)*W1(1) do j=2,nlobit cm(1)=cm(1)+z(2,j)*W1(j) cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1))) enddo cm(1)=cm(1)/sumW(nlobit) cm(2)=cm(2)/sumW(nlobit) if (cm(1).gt.Big .or. cm(1).lt.-Big .or. & cm(2).gt.Big .or. cm(2).lt.-Big) then cd write (iout,'(a)') cd & 'Unexpected error in GenSide - CM coordinates too large.' cd write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2) cd write (*,'(a)') cd & 'Unexpected error in GenSide - CM coordinates too large.' cd write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2) fail=.true. return endif cd print *,'CM:',cm(1),cm(2) C C Find the largest search distance from CM C radmax=0.0D0 do i=1,nlobit do j=2,3 do k=2,3 a(j-1,k-1)=gaussc(j,k,i,it) enddo enddo #ifdef NAG call f02faf('N','U',2,a,3,eig,work,100,ifail) #else call djacob(2,3,10000,1.0d-10,a,vec,eig) #endif #ifdef MPI if (lprint) then print *,'*************** CG Processor',me print *,'CM:',cm(1),cm(2) write (iout,*) '*************** CG Processor',me write (iout,*) 'CM:',cm(1),cm(2) print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2) write (iout,'(A,8f10.5)') & 'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2) endif #endif if (eig(1).lt.eig_limit) then write(iout,'(a)') & 'From Mult_Norm: Eigenvalues of A are too small.' write(*,'(a)') & 'From Mult_Norm: Eigenvalues of A are too small.' fail=.true. return endif radius=0.0D0 cd print *,'i=',i do j=1,2 radius=radius+pinorm(z(j+1,i)-cm(j))**2 enddo radius=dsqrt(radius)+3.0D0/dsqrt(eig(1)) if (radius.gt.radmax) radmax=radius enddo if (radmax.gt.pi) radmax=pi C C Determine the boundaries of the search rectangle. C if (lprint) then print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) ) print '(a,4(1pe14.4))','radmax: ',radmax endif box(1,1)=dmax1(cm(1)-radmax,0.0D0) box(2,1)=dmin1(cm(1)+radmax,pi) box(1,2)=cm(2)-radmax box(2,2)=cm(2)+radmax if (lprint) then #ifdef MPI print *,'CG Processor',me,' Array BOX:' #else print *,'Array BOX:' #endif print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2) print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) ) #ifdef MPI write (iout,*)'CG Processor',me,' Array BOX:' #else write (iout,*)'Array BOX:' #endif write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2) write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) ) endif if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then #ifdef MPI write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.' write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.' #else c write (iout,'(a)') 'Bad sampling box.' #endif fail=.true. return endif which_lobe=ran_number(0.0D0,sumW(nlobit)) c print '(a,1pe14.4)','which_lobe=',which_lobe do i=1,nlobit if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1 enddo 1 ilob=i c print *,'ilob=',ilob,' nlob=',nlob(it) do i=2,3 cm(i-1)=z(i,ilob) do j=2,3 a(i-1,j-1)=gaussc(i,j,ilob,it) enddo enddo cd print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.' call mult_norm1(3,2,a,cm,box,y,fail) if (fail) return al=y(1) om=pinorm(y(2)) cd print *,'al=',al,' om=',om cd stop return end c--------------------------------------------------------------------------- double precision function ran_number(x1,x2) C Calculate a random real number from the range (x1,x2). implicit real*8 (a-h,o-z) include 'DIMENSIONS' double precision x1,x2,fctor data fctor /2147483647.0D0/ #ifdef MPI include "mpif.h" include 'COMMON.SETUP' ran_number=x1+(x2-x1)*prng_next(me) #else call vrnd(ix,1) ran_number=x1+(x2-x1)*ix/fctor #endif return end c-------------------------------------------------------------------------- integer function iran_num(n1,n2) C Calculate a random integer number from the range (n1,n2). implicit real*8 (a-h,o-z) include 'DIMENSIONS' integer n1,n2,ix real fctor /2147483647.0/ #ifdef MPI include "mpif.h" include 'COMMON.SETUP' ix=n1+(n2-n1+1)*prng_next(me) if (ix.lt.n1) ix=n1 if (ix.gt.n2) ix=n2 iran_num=ix #else call vrnd(ix,1) ix=n1+(n2-n1+1)*(ix/fctor) if (ix.gt.n2) ix=n2 iran_num=ix #endif return end c-------------------------------------------------------------------------- double precision function binorm(x1,x2,sigma1,sigma2,ak) implicit real*8 (a-h,o-z) c print '(a)','Enter BINORM.' alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2) aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2) seg=sigma1/(sigma1+ak*sigma2) alen=ran_number(0.0D0,1.0D0) if (alen.lt.seg) then binorm=anorm_distr(x1,sigma1,alowb,aupb) else binorm=anorm_distr(x2,sigma2,alowb,aupb) endif c print '(a)','Exiting BINORM.' return end c----------------------------------------------------------------------- c double precision function anorm_distr(x,sigma,alowb,aupb) c implicit real*8 (a-h,o-z) c print '(a)','Enter ANORM_DISTR.' c 10 y=ran_number(alowb,aupb) c expon=dexp(-0.5D0*((y-x)/sigma)**2) c ran=ran_number(0.0D0,1.0D0) c if (expon.lt.ran) goto 10 c anorm_distr=y c print '(a)','Exiting ANORM_DISTR.' c return c end c----------------------------------------------------------------------- double precision function anorm_distr(x,sigma,alowb,aupb) implicit real*8 (a-h,o-z) c to make a normally distributed deviate with zero mean and unit variance c integer iset real fac,gset,rsq,v1,v2,ran1 save iset,gset data iset/0/ if(iset.eq.0) then 1 v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0 rsq=v1**2+v2**2 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1 fac=sqrt(-2.0d0*log(rsq)/rsq) gset=v1*fac gaussdev=v2*fac iset=1 else gaussdev=gset iset=0 endif anorm_distr=x+gaussdev*sigma return end c------------------------------------------------------------------------ subroutine mult_norm(lda,n,a,x,fail) C C Generate the vector X whose elements obey the multiple-normal distribution C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A, C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest C eigenvalue of the matrix A is close to 0. C implicit double precision (a-h,o-z) double precision a(lda,n),x(n),eig(100),vec(3,3),work(100) double precision eig_limit /1.0D-8/ logical fail fail=.false. c print '(a)','Enter MULT_NORM.' C C Find the smallest eigenvalue of the matrix A. C c do i=1,n c print '(8f10.5)',(a(i,j),j=1,n) c enddo #ifdef NAG call f02faf('V','U',2,a,lda,eig,work,100,ifail) #else call djacob(2,lda,10000,1.0d-10,a,vec,eig) #endif c print '(8f10.5)',(eig(i),i=1,n) C print '(a)' c do i=1,n c print '(8f10.5)',(a(i,j),j=1,n) c enddo if (eig(1).lt.eig_limit) then print *,'From Mult_Norm: Eigenvalues of A are too small.' fail=.true. return endif C C Generate points following the normal distributions along the principal C axes of the moment matrix. Store in WORK. C do i=1,n sigma=1.0D0/dsqrt(eig(i)) alim=-3.0D0*sigma work(i)=anorm_distr(0.0D0,sigma,-alim,alim) enddo C C Transform the vector of normal variables back to the original basis. C do i=1,n xi=0.0D0 do j=1,n xi=xi+a(i,j)*work(j) enddo x(i)=xi enddo return end c------------------------------------------------------------------------ subroutine mult_norm1(lda,n,a,z,box,x,fail) C C Generate the vector X whose elements obey the multi-gaussian multi-dimensional C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the C leading dimension of the moment matrix A, n is the dimension of the C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the C smallest eigenvalue of the matrix A is close to 0. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif double precision a(lda,n),z(n),x(n),box(n,n) double precision etmp include 'COMMON.IOUNITS' #ifdef MP include 'COMMON.SETUP' #endif logical fail C C Generate points following the normal distributions along the principal C axes of the moment matrix. Store in WORK. C cd print *,'CG Processor',me,' entered MultNorm1.' cd print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2) cd do i=1,n cd print *,i,box(1,i),box(2,i) cd enddo istep = 0 10 istep = istep + 1 if (istep.gt.10000) then c write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps', c & ' in MultNorm1.' c write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps', c & ' in MultNorm1.' c write (iout,*) 'box',box c write (iout,*) 'a',a c write (iout,*) 'z',z fail=.true. return endif do i=1,n x(i)=ran_number(box(1,i),box(2,i)) enddo ww=0.0D0 do i=1,n xi=pinorm(x(i)-z(i)) ww=ww+0.5D0*a(i,i)*xi*xi do j=i+1,n ww=ww+a(i,j)*xi*pinorm(x(j)-z(j)) enddo enddo dec=ran_number(0.0D0,1.0D0) c print *,(x(i),i=1,n),ww,dexp(-ww),dec crc if (dec.gt.dexp(-ww)) goto 10 if(-ww.lt.100) then etmp=dexp(-ww) else return endif if (dec.gt.etmp) goto 10 cd print *,'CG Processor',me,' exitting MultNorm1.' return end c crc-------------------------------------- subroutine overlap_sc(scfail) c Internal and cartesian coordinates must be consistent as input, c and will be up-to-date on return. c At the end of this procedure, scfail is true if there are c overlapping residues left, or false otherwise (success) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.VAR' include 'COMMON.SBRIDGE' include 'COMMON.IOUNITS' logical had_overlaps,fail,scfail integer ioverlap(maxres),ioverlap_last had_overlaps=.false. call overlap_sc_list(ioverlap,ioverlap_last) if (ioverlap_last.gt.0) then write (iout,*) '#OVERLAPing residues ',ioverlap_last write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last) had_overlaps=.true. endif maxsi=1000 do k=1,1000 if (ioverlap_last.eq.0) exit do ires=1,ioverlap_last i=ioverlap(ires) iti=iabs(itype(i)) if (iti.ne.10 .and. iti.lt.ntyp1 .and.iti.gt.0) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) call gen_side(iti,theta(i+1),alph(i),omeg(i),fail) nsi=nsi+1 enddo if(fail) goto 999 endif enddo call chainbuild call overlap_sc_list(ioverlap,ioverlap_last) c write (iout,*) 'Overlaping residues ',ioverlap_last, c & (ioverlap(j),j=1,ioverlap_last) enddo if (k.le.1000.and.ioverlap_last.eq.0) then scfail=.false. if (had_overlaps) then write (iout,*) '#OVERLAPing all corrected after ',k, & ' random generation' endif else scfail=.true. write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last) endif return 999 continue write (iout,'(a30,i5,a12,i4)') & '#OVERLAP FAIL in gen_side after',maxsi, & 'iter for RES',i scfail=.true. return end subroutine overlap_sc_list(ioverlap,ioverlap_last) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.IOUNITS' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.VAR' include 'COMMON.CALC' logical fail integer ioverlap(maxres),ioverlap_last data redfac /0.5D0/ ioverlap_last=0 C Check for SC-SC overlaps and mark residues c print *,'>>overlap_sc nnt=',nnt,' nct=',nct ind=0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) itypi1=iabs(itype(i+1)) if (itypi.gt.ntyp .or. itypi.le.0) cycle xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) dsci_inv=dsc_inv(itypi) c do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.gt.ntyp .or. itypj.le.0) cycle dscj_inv=dsc_inv(itypj) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) chi12=chi1*chi2 chip1=chip(itypi) chip2=chip(itypj) chip12=chip1*chip2 alf1=alp(itypi) alf2=alp(itypj) alf12=0.5D0*(alf1+alf2) if (j.gt.i+1) then rcomp=sigmaii(itypi,itypj) else rcomp=sigma(itypi,itypj) endif c print '(2(a3,2i3),a3,2f10.5)', c & ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j) c & ,rcomp xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) call sc_angular sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) rij_shift=1.0D0/rij-sig+sig0ij ct if ( 1.0/rij .lt. redfac*rcomp .or. ct & rij_shift.le.0.0D0 ) then if ( rij_shift.le.0.0D0 ) then cd write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)') cd & 'overlap SC-SC: i=',i,' j=',j, cd & ' dist=',dist(nres+i,nres+j),' rcomp=', cd & rcomp,1.0/rij,rij_shift ioverlap_last=ioverlap_last+1 ioverlap(ioverlap_last)=i do k=1,ioverlap_last-1 if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1 enddo ioverlap_last=ioverlap_last+1 ioverlap(ioverlap_last)=j do k=1,ioverlap_last-1 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1 enddo endif enddo enddo enddo return end