X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_CSA_DiL%2Fgen_rand_conf.F;fp=source%2Funres%2Fsrc_CSA_DiL%2Fgen_rand_conf.F;h=0000000000000000000000000000000000000000;hp=78d4ccae7af9b2057cca2af74651f401e98f93eb;hb=de4bc5453ea46e111d936cb85e1758ed21c08fcd;hpb=b75425747e3e2b448ca5e0ef8367712e1f339124 diff --git a/source/unres/src_CSA_DiL/gen_rand_conf.F b/source/unres/src_CSA_DiL/gen_rand_conf.F deleted file mode 100644 index 78d4cca..0000000 --- a/source/unres/src_CSA_DiL/gen_rand_conf.F +++ /dev/null @@ -1,911 +0,0 @@ - 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)),abs(itype(3))) -c write(iout,*)'phi(4)=',rad2deg*phi(4) - if (nstart.lt.3) theta(3)=gen_theta(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=abs(itype(i-1)) - it2=abs(itype(i-2)) - it=abs(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) - 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=abs(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=abs(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=abs(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=abs(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 - 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=abs(itype(i)) - itypi1=abs(itype(i+1)) - 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=itype(j) - 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