--- /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=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
--- /dev/null
+ subroutine sc_move(n_start,n_end,n_maxtry,e_drop,
+ + n_fun,etot)
+c Perform a quick search over side-chain arrangments (over
+c residues n_start to n_end) for a given (frozen) CA trace
+c Only side-chains are minimized (at most n_maxtry times each),
+c not CA positions
+c Stops if energy drops by e_drop, otherwise tries all residues
+c in the given range
+c If there is an energy drop, full minimization may be useful
+c n_start, n_end CAN be modified by this routine, but only if
+c out of bounds (n_start <= 1, n_end >= nres, n_start < n_end)
+c NOTE: this move should never increase the energy
+crc implicit none
+
+c Includes
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.HEADER'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.FFIELD'
+
+c External functions
+ integer iran_num
+ external iran_num
+
+c Input arguments
+ integer n_start,n_end,n_maxtry
+ double precision e_drop
+
+c Output arguments
+ integer n_fun
+ double precision etot
+
+c Local variables
+ double precision energy(0:n_ene)
+ double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
+ double precision orig_e,cur_e
+ integer n,n_steps,n_first,n_cur,n_tot,i
+ double precision orig_w(n_ene)
+ double precision wtime
+
+
+c Set non side-chain weights to zero (minimization is faster)
+c NOTE: e(2) does not actually depend on the side-chain, only CA
+ orig_w(2)=wscp
+ orig_w(3)=welec
+ orig_w(4)=wcorr
+ orig_w(5)=wcorr5
+ orig_w(6)=wcorr6
+ orig_w(7)=wel_loc
+ orig_w(8)=wturn3
+ orig_w(9)=wturn4
+ orig_w(10)=wturn6
+ orig_w(11)=wang
+ orig_w(13)=wtor
+ orig_w(14)=wtor_d
+ orig_w(15)=wvdwpp
+
+ wscp=0.D0
+ welec=0.D0
+ wcorr=0.D0
+ wcorr5=0.D0
+ wcorr6=0.D0
+ wel_loc=0.D0
+ wturn3=0.D0
+ wturn4=0.D0
+ wturn6=0.D0
+ wang=0.D0
+ wtor=0.D0
+ wtor_d=0.D0
+ wvdwpp=0.D0
+
+c Make sure n_start, n_end are within proper range
+ if (n_start.lt.2) n_start=2
+ if (n_end.gt.nres-1) n_end=nres-1
+crc if (n_start.lt.n_end) then
+ if (n_start.gt.n_end) then
+ n_start=2
+ n_end=nres-1
+ endif
+
+c Save the initial values of energy and coordinates
+cd call chainbuild
+cd call etotal(energy)
+cd write (iout,*) 'start sc ene',energy(0)
+cd call enerprint(energy(0))
+crc etot=energy(0)
+ n_fun=0
+crc orig_e=etot
+crc cur_e=orig_e
+crc do i=2,nres-1
+crc cur_alph(i)=alph(i)
+crc cur_omeg(i)=omeg(i)
+crc enddo
+
+ct wtime=MPI_WTIME()
+c Try (one by one) all specified residues, starting from a
+c random position in sequence
+c Stop early if the energy has decreased by at least e_drop
+ n_tot=n_end-n_start+1
+ n_first=iran_num(0,n_tot-1)
+ n_steps=0
+ n=0
+crc do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop)
+ do while (n.lt.n_tot)
+ n_cur=n_start+mod(n_first+n,n_tot)
+ call single_sc_move(n_cur,n_maxtry,e_drop,
+ + n_steps,n_fun,etot)
+c If a lower energy was found, update the current structure...
+crc if (etot.lt.cur_e) then
+crc cur_e=etot
+crc do i=2,nres-1
+crc cur_alph(i)=alph(i)
+crc cur_omeg(i)=omeg(i)
+crc enddo
+crc else
+c ...else revert to the previous one
+crc etot=cur_e
+crc do i=2,nres-1
+crc alph(i)=cur_alph(i)
+crc omeg(i)=cur_omeg(i)
+crc enddo
+crc endif
+ n=n+1
+cd
+cd call chainbuild
+cd call etotal(energy)
+cd print *,'running',n,energy(0)
+ enddo
+
+cd call chainbuild
+cd call etotal(energy)
+cd write (iout,*) 'end sc ene',energy(0)
+
+c Put the original weights back to calculate the full energy
+ wscp=orig_w(2)
+ welec=orig_w(3)
+ wcorr=orig_w(4)
+ wcorr5=orig_w(5)
+ wcorr6=orig_w(6)
+ wel_loc=orig_w(7)
+ wturn3=orig_w(8)
+ wturn4=orig_w(9)
+ wturn6=orig_w(10)
+ wang=orig_w(11)
+ wtor=orig_w(13)
+ wtor_d=orig_w(14)
+ wvdwpp=orig_w(15)
+
+crc n_fun=n_fun+1
+ct write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
+ return
+ end
+
+c-------------------------------------------------------------
+
+ subroutine single_sc_move(res_pick,n_maxtry,e_drop,
+ + n_steps,n_fun,e_sc)
+c Perturb one side-chain (res_pick) and minimize the
+c neighbouring region, keeping all CA's and non-neighbouring
+c side-chains fixed
+c Try until e_drop energy improvement is achieved, or n_maxtry
+c attempts have been made
+c At the start, e_sc should contain the side-chain-only energy(0)
+c nsteps and nfun for this move are ADDED to n_steps and n_fun
+crc implicit none
+
+c Includes
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.CHAIN'
+ include 'COMMON.MINIM'
+ include 'COMMON.FFIELD'
+ include 'COMMON.IOUNITS'
+
+c External functions
+ double precision dist
+ external dist
+
+c Input arguments
+ integer res_pick,n_maxtry
+ double precision e_drop
+
+c Input/Output arguments
+ integer n_steps,n_fun
+ double precision e_sc
+
+c Local variables
+ logical fail
+ integer i,j
+ integer nres_moved
+ integer iretcode,loc_nfun,orig_maxfun,n_try
+ double precision sc_dist,sc_dist_cutoff
+ double precision energy(0:n_ene),orig_e,cur_e
+ double precision evdw,escloc
+ double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
+ double precision var(maxvar)
+
+ double precision orig_theta(1:nres),orig_phi(1:nres),
+ + orig_alph(1:nres),orig_omeg(1:nres)
+
+
+c Define what is meant by "neighbouring side-chain"
+ sc_dist_cutoff=5.0D0
+
+c Don't do glycine or ends
+ i=itype(res_pick)
+ if (i.eq.10 .or. i.eq.21) return
+
+c Freeze everything (later will relax only selected side-chains)
+ mask_r=.true.
+ do i=1,nres
+ mask_phi(i)=0
+ mask_theta(i)=0
+ mask_side(i)=0
+ enddo
+
+c Find the neighbours of the side-chain to move
+c and save initial variables
+crc orig_e=e_sc
+crc cur_e=orig_e
+ nres_moved=0
+ do i=2,nres-1
+c Don't do glycine (itype(j)==10)
+ if (itype(i).ne.10) then
+ sc_dist=dist(nres+i,nres+res_pick)
+ else
+ sc_dist=sc_dist_cutoff
+ endif
+ if (sc_dist.lt.sc_dist_cutoff) then
+ nres_moved=nres_moved+1
+ mask_side(i)=1
+ cur_alph(i)=alph(i)
+ cur_omeg(i)=omeg(i)
+ endif
+ enddo
+
+ call chainbuild
+ call egb1(evdw)
+ call esc(escloc)
+ e_sc=wsc*evdw+wscloc*escloc
+cd call etotal(energy)
+cd print *,'new ',(energy(k),k=0,n_ene)
+ orig_e=e_sc
+ cur_e=orig_e
+
+ n_try=0
+ do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
+c Move the selected residue (don't worry if it fails)
+ call gen_side(itype(res_pick),theta(res_pick+1),
+ + alph(res_pick),omeg(res_pick),fail)
+
+c Minimize the side-chains starting from the new arrangement
+ call geom_to_var(nvar,var)
+ orig_maxfun=maxfun
+ maxfun=7
+
+crc do i=1,nres
+crc orig_theta(i)=theta(i)
+crc orig_phi(i)=phi(i)
+crc orig_alph(i)=alph(i)
+crc orig_omeg(i)=omeg(i)
+crc enddo
+
+ call minimize_sc1(e_sc,var,iretcode,loc_nfun)
+
+cv write(*,'(2i3,2f12.5,2i3)')
+cv & res_pick,nres_moved,orig_e,e_sc-cur_e,
+cv & iretcode,loc_nfun
+
+c$$$ if (iretcode.eq.8) then
+c$$$ write(iout,*)'Coordinates just after code 8'
+c$$$ call chainbuild
+c$$$ call all_varout
+c$$$ call flush(iout)
+c$$$ do i=1,nres
+c$$$ theta(i)=orig_theta(i)
+c$$$ phi(i)=orig_phi(i)
+c$$$ alph(i)=orig_alph(i)
+c$$$ omeg(i)=orig_omeg(i)
+c$$$ enddo
+c$$$ write(iout,*)'Coordinates just before code 8'
+c$$$ call chainbuild
+c$$$ call all_varout
+c$$$ call flush(iout)
+c$$$ endif
+
+ n_fun=n_fun+loc_nfun
+ maxfun=orig_maxfun
+ call var_to_geom(nvar,var)
+
+c If a lower energy was found, update the current structure...
+ if (e_sc.lt.cur_e) then
+cv call chainbuild
+cv call etotal(energy)
+cd call egb1(evdw)
+cd call esc(escloc)
+cd e_sc1=wsc*evdw+wscloc*escloc
+cd print *,' new',e_sc1,energy(0)
+cv print *,'new ',energy(0)
+cd call enerprint(energy(0))
+ cur_e=e_sc
+ do i=2,nres-1
+ if (mask_side(i).eq.1) then
+ cur_alph(i)=alph(i)
+ cur_omeg(i)=omeg(i)
+ endif
+ enddo
+ else
+c ...else revert to the previous one
+ e_sc=cur_e
+ do i=2,nres-1
+ if (mask_side(i).eq.1) then
+ alph(i)=cur_alph(i)
+ omeg(i)=cur_omeg(i)
+ endif
+ enddo
+ endif
+ n_try=n_try+1
+
+ enddo
+ n_steps=n_steps+n_try
+
+c Reset the minimization mask_r to false
+ mask_r=.false.
+
+ return
+ end
+
+c-------------------------------------------------------------
+
+ subroutine sc_minimize(etot,iretcode,nfun)
+c Minimizes side-chains only, leaving backbone frozen
+crc implicit none
+
+c Includes
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.FFIELD'
+
+c Output arguments
+ double precision etot
+ integer iretcode,nfun
+
+c Local variables
+ integer i
+ double precision orig_w(n_ene),energy(0:n_ene)
+ double precision var(maxvar)
+
+
+c Set non side-chain weights to zero (minimization is faster)
+c NOTE: e(2) does not actually depend on the side-chain, only CA
+ orig_w(2)=wscp
+ orig_w(3)=welec
+ orig_w(4)=wcorr
+ orig_w(5)=wcorr5
+ orig_w(6)=wcorr6
+ orig_w(7)=wel_loc
+ orig_w(8)=wturn3
+ orig_w(9)=wturn4
+ orig_w(10)=wturn6
+ orig_w(11)=wang
+ orig_w(13)=wtor
+ orig_w(14)=wtor_d
+
+ wscp=0.D0
+ welec=0.D0
+ wcorr=0.D0
+ wcorr5=0.D0
+ wcorr6=0.D0
+ wel_loc=0.D0
+ wturn3=0.D0
+ wturn4=0.D0
+ wturn6=0.D0
+ wang=0.D0
+ wtor=0.D0
+ wtor_d=0.D0
+
+c Prepare to freeze backbone
+ do i=1,nres
+ mask_phi(i)=0
+ mask_theta(i)=0
+ mask_side(i)=1
+ enddo
+
+c Minimize the side-chains
+ mask_r=.true.
+ call geom_to_var(nvar,var)
+ call minimize(etot,var,iretcode,nfun)
+ call var_to_geom(nvar,var)
+ mask_r=.false.
+
+c Put the original weights back and calculate the full energy
+ wscp=orig_w(2)
+ welec=orig_w(3)
+ wcorr=orig_w(4)
+ wcorr5=orig_w(5)
+ wcorr6=orig_w(6)
+ wel_loc=orig_w(7)
+ wturn3=orig_w(8)
+ wturn4=orig_w(9)
+ wturn6=orig_w(10)
+ wang=orig_w(11)
+ wtor=orig_w(13)
+ wtor_d=orig_w(14)
+
+ call chainbuild
+ call etotal(energy)
+ etot=energy(0)
+
+ return
+ end
+
+c-------------------------------------------------------------
+ subroutine minimize_sc1(etot,x,iretcode,nfun)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
+ include 'COMMON.IOUNITS'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.MINIM'
+ common /srutu/ icall
+ dimension iv(liv)
+ double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
+ double precision energia(0:n_ene)
+ external func,gradient,fdum
+ external func_restr1,grad_restr1
+ logical not_done,change,reduce
+ common /przechowalnia/ v
+
+ call deflt(2,iv,liv,lv,v)
+* 12 means fresh start, dont call deflt
+ iv(1)=12
+* max num of fun calls
+ if (maxfun.eq.0) maxfun=500
+ iv(17)=maxfun
+* max num of iterations
+ if (maxmin.eq.0) maxmin=1000
+ iv(18)=maxmin
+* controls output
+ iv(19)=2
+* selects output unit
+c iv(21)=iout
+ iv(21)=0
+* 1 means to print out result
+ iv(22)=0
+* 1 means to print out summary stats
+ iv(23)=0
+* 1 means to print initial x and d
+ iv(24)=0
+* min val for v(radfac) default is 0.1
+ v(24)=0.1D0
+* max val for v(radfac) default is 4.0
+ v(25)=2.0D0
+c v(25)=4.0D0
+* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
+* the sumsl default is 0.1
+ v(26)=0.1D0
+* false conv if (act fnctn decrease) .lt. v(34)
+* the sumsl default is 100*machep
+ v(34)=v(34)/100.0D0
+* absolute convergence
+ if (tolf.eq.0.0D0) tolf=1.0D-4
+ v(31)=tolf
+* relative convergence
+ if (rtolf.eq.0.0D0) rtolf=1.0D-4
+ v(32)=rtolf
+* controls initial step size
+ v(35)=1.0D-1
+* large vals of d correspond to small components of step
+ do i=1,nphi
+ d(i)=1.0D-1
+ enddo
+ do i=nphi+1,nvar
+ d(i)=1.0D-1
+ enddo
+ IF (mask_r) THEN
+ call x2xx(x,xx,nvar_restr)
+ call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,
+ & iv,liv,lv,v,idum,rdum,fdum)
+ call xx2x(x,xx)
+ ELSE
+ call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
+ ENDIF
+ etot=v(10)
+ iretcode=iv(1)
+ nfun=iv(6)
+
+ return
+ end
+************************************************************************
+ subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.FFIELD'
+ include 'COMMON.INTERACT'
+ include 'COMMON.TIME1'
+ common /chuju/ jjj
+ double precision energia(0:n_ene),evdw,escloc
+ integer jjj
+ double precision ufparm,e1,e2
+ external ufparm
+ integer uiparm(1)
+ real*8 urparm(1)
+ dimension x(maxvar)
+ nfl=nf
+ icg=mod(nf,2)+1
+
+#ifdef OSF
+c Intercept NaNs in the coordinates, before calling etotal
+ x_sum=0.D0
+ do i=1,n
+ x_sum=x_sum+x(i)
+ enddo
+ FOUND_NAN=.false.
+ if (x_sum.ne.x_sum) then
+ write(iout,*)" *** func_restr1 : Found NaN in coordinates"
+ f=1.0D+73
+ FOUND_NAN=.true.
+ return
+ endif
+#endif
+
+ call var_to_geom_restr(n,x)
+ call zerograd
+ call chainbuild
+cd write (iout,*) 'ETOTAL called from FUNC'
+ call egb1(evdw)
+ call esc(escloc)
+ f=wsc*evdw+wscloc*escloc
+cd call etotal(energia(0))
+cd f=wsc*energia(1)+wscloc*energia(12)
+cd print *,f,evdw,escloc,energia(0)
+C
+C Sum up the components of the Cartesian gradient.
+C
+ do i=1,nct
+ do j=1,3
+ gradx(j,i,icg)=wsc*gvdwx(j,i)
+ enddo
+ enddo
+
+ return
+ end
+c-------------------------------------------------------
+ subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.FFIELD'
+ include 'COMMON.IOUNITS'
+ external ufparm
+ integer uiparm(1)
+ double precision urparm(1)
+ dimension x(maxvar),g(maxvar)
+
+ icg=mod(nf,2)+1
+ if (nf-nfl+1) 20,30,40
+ 20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
+c write (iout,*) 'grad 20'
+ if (nf.eq.0) return
+ goto 40
+ 30 call var_to_geom_restr(n,x)
+ call chainbuild
+C
+C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
+C
+ 40 call cartder
+C
+C Convert the Cartesian gradient into internal-coordinate gradient.
+C
+
+ ig=0
+ ind=nres-2
+ do i=2,nres-2
+ IF (mask_phi(i+2).eq.1) THEN
+ gphii=0.0D0
+ do j=i+1,nres-1
+ ind=ind+1
+ do k=1,3
+ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
+ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
+ enddo
+ enddo
+ ig=ig+1
+ g(ig)=gphii
+ ELSE
+ ind=ind+nres-1-i
+ ENDIF
+ enddo
+
+
+ ind=0
+ do i=1,nres-2
+ IF (mask_theta(i+2).eq.1) THEN
+ ig=ig+1
+ gthetai=0.0D0
+ do j=i+1,nres-1
+ ind=ind+1
+ do k=1,3
+ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
+ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
+ enddo
+ enddo
+ g(ig)=gthetai
+ ELSE
+ ind=ind+nres-1-i
+ ENDIF
+ enddo
+
+ do i=2,nres-1
+ if (itype(i).ne.10) then
+ IF (mask_side(i).eq.1) THEN
+ ig=ig+1
+ galphai=0.0D0
+ do k=1,3
+ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
+ enddo
+ g(ig)=galphai
+ ENDIF
+ endif
+ enddo
+
+
+ do i=2,nres-1
+ if (itype(i).ne.10) then
+ IF (mask_side(i).eq.1) THEN
+ ig=ig+1
+ gomegai=0.0D0
+ do k=1,3
+ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
+ enddo
+ g(ig)=gomegai
+ ENDIF
+ endif
+ enddo
+
+C
+C Add the components corresponding to local energy terms.
+C
+
+ ig=0
+ igall=0
+ do i=4,nres
+ igall=igall+1
+ if (mask_phi(i).eq.1) then
+ ig=ig+1
+ g(ig)=g(ig)+gloc(igall,icg)
+ endif
+ enddo
+
+ do i=3,nres
+ igall=igall+1
+ if (mask_theta(i).eq.1) then
+ ig=ig+1
+ g(ig)=g(ig)+gloc(igall,icg)
+ endif
+ enddo
+
+ do ij=1,2
+ do i=2,nres-1
+ if (itype(i).ne.10) then
+ igall=igall+1
+ if (mask_side(i).eq.1) then
+ ig=ig+1
+ g(ig)=g(ig)+gloc(igall,icg)
+ endif
+ endif
+ enddo
+ enddo
+
+cd do i=1,ig
+cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
+cd enddo
+ return
+ end
+C-----------------------------------------------------------------------------
+ subroutine egb1(evdw)
+C
+C This subroutine calculates the interaction energy of nonbonded side chains
+C assuming the Gay-Berne potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ logical lprn
+ evdw=0.0D0
+c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ lprn=.false.
+c if (icall.eq.0) lprn=.true.
+ 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
+C Calculate SC interaction energy.
+C
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
+ 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)
+C For diagnostics only!!!
+c chi1=0.0D0
+c chi2=0.0D0
+c chi12=0.0D0
+c chip1=0.0D0
+c chip2=0.0D0
+c chip12=0.0D0
+c alf1=0.0D0
+c alf2=0.0D0
+c alf12=0.0D0
+ 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)
+C Calculate angle-dependent terms of energy and contributions to their
+C derivatives.
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+C I hate to put IF's in the loops, but here don't have another choice!!!!
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+cd & restyp(itypi),i,restyp(itypj),j,
+cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
+ return
+ endif
+ sigder=-sig*sigsq
+c---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+cd & restyp(itypi),i,restyp(itypj),j,
+cd & epsi,sigm,chi1,chi2,chip1,chip2,
+cd & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+cd & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+cd & evdwij
+ endif
+
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
+ & 'evdw',i,j,evdwij
+
+C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+C Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+C Calculate angular part of the gradient.
+ call sc_grad
+ ENDIF
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ end
+C-----------------------------------------------------------------------------