+++ /dev/null
- 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
projects / unres.git / blobdiff