subroutine gen_rand_conf(nstart,*) 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,back,fail integer nstart integer i,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit double precision gen_theta,gen_phi,dist 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 .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 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 .and. it1.ne.ntyp1) 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 none include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.FFIELD' double precision redfac /0.5D0/ integer i,j,k,iti,itj,iteli,itelj double precision rcomp double precision dist 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.TORCNSTR" include 'COMMON.GEO' include 'COMMON.BOUNDS' if (raw_psipred .or. ndih_constr.eq.0) then gen_phi=ran_number(-pi,pi) else C 8/13/98 Generate phi using pre-defined boundaries gen_phi=ran_number(phibound(1,i),phibound(2,i)) endif 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,lprn /.false./ 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 if (lprn) then write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.' write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.' endif #else if (lprn) 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,'(18i5)') (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) 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_extconf call overlap_sc_list(ioverlap,ioverlap_last) write (iout,*) 'Overlaping residues ',ioverlap_last, & (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/ write (iout,*) "overlap_sc_list" c write(iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e write(iout,*) "nnt",nnt," nct",nct ioverlap_last=0 C Check for SC-SC overlaps and mark residues c print *,'>>overlap_sc nnt=',nnt,' nct=',nct ind=0 c do i=iatsc_s,iatsc_e do i=nnt,nct itypi=iabs(itype(i)) itypi1=iabs(itype(i+1)) if (itypi.eq.ntyp1) 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 c do iint=1,nint_gr(i) c do j=istart(i,iint),iend(i,iint) do j=i+1,nct ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle c write (iout,*) "i,j",i,j," itypi,itypj",itypi,itypj 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 write (iout,'(2(a3,2i5),a3,2f10.5)'), c & ' i=',i,itypi,' j=',j,itypj,' 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 c write (iout,*) "dxj",dxj," dyj",dyj," dzj",dzj c write (iout,*) "erij",erij c write (iout,*) "om1",om1," om2",om2," om12",om12, c & " faceps1",faceps1," eps1",eps1 c write (iout,*) "sigsq",sigsq sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) rij_shift=1.0D0/rij-sig+sig0ij c write (iout,*) "rij",1.0d0/rij," rij_shift",rij_shift, c & " sig",sig," sig0ij",sig0ij c if ( rij_shift.le.0.0D0 ) then if ( rij_shift/sig0ij.le.0.1D0 ) then c write (iout,*) "overlap",i,j write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)') & 'overlap SC-SC: i=',i,' j=',j, & ' dist=',dist(nres+i,nres+j),' rcomp=', & 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 c enddo enddo return end