X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_Eshel%2Fgen_rand_conf.F;fp=source%2Funres%2Fsrc_Eshel%2Fgen_rand_conf.F;h=6cc31ba6e4c2cc5e264bccb32e90a9d575afb602;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_Eshel/gen_rand_conf.F b/source/unres/src_Eshel/gen_rand_conf.F new file mode 100644 index 0000000..6cc31ba --- /dev/null +++ b/source/unres/src_Eshel/gen_rand_conf.F @@ -0,0 +1,910 @@ + 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=itype(2) + phi(4)=gen_phi(4,itype(2),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=itype(i-1) + it2=itype(i-2) + it=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=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=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=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)*y(k)+bthet(k,it)*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=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=itype(i) + itypi1=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