source/unres/src_CSA_DiL/gen_rand_conf.F

   1       subroutine gen_rand_conf(nstart,*)
   2 C Generate random conformation or chain cut and regrowth.
   3       implicit real*8 (a-h,o-z)
   4       include 'DIMENSIONS'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.LOCAL'
   7       include 'COMMON.VAR'
   8       include 'COMMON.INTERACT'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.MCM'
  11       include 'COMMON.GEO'
  12       include 'COMMON.CONTROL'
  13       logical overlap,back,fail
  14 cd    print *,' CG Processor',me,' maxgen=',maxgen
  15       maxsi=100
  16 cd    write (iout,*) 'Gen_Rand_conf: nstart=',nstart
  17       if (nstart.lt.5) then
  18         it1=iabs(itype(2))
  19         phi(4)=gen_phi(4,iabs(itype(2)),abs(itype(3)))
  20 c       write(iout,*)'phi(4)=',rad2deg*phi(4)
  21         if (nstart.lt.3) theta(3)=gen_theta(itype(2),pi,phi(4))
  22 c       write(iout,*)'theta(3)=',rad2deg*theta(3)
  23         if (it1.ne.10) then
  24           nsi=0
  25           fail=.true.
  26           do while (fail.and.nsi.le.maxsi)
  27             call gen_side(it1,theta(3),alph(2),omeg(2),fail)
  28             nsi=nsi+1
  29           enddo
  30           if (nsi.gt.maxsi) return1
  31         endif ! it1.ne.10
  32         call orig_frame
  33         i=4
  34         nstart=4
  35       else
  36         i=nstart
  37         nstart=max0(i,4)
  38       endif
  39
  40       maxnit=0
  41
  42       nit=0
  43       niter=0
  44       back=.false.
  45       do while (i.le.nres .and. niter.lt.maxgen)
  46         if (i.lt.nstart) then
  47           if(iprint.gt.1) then
  48           write (iout,'(/80(1h*)/2a/80(1h*))')
  49      &          'Generation procedure went down to ',
  50      &          'chain beginning. Cannot continue...'
  51           write (*,'(/80(1h*)/2a/80(1h*))')
  52      &          'Generation procedure went down to ',
  53      &          'chain beginning. Cannot continue...'
  54           endif
  55           return1
  56         endif
  57         it1=abs(itype(i-1))
  58         it2=abs(itype(i-2))
  59         it=abs(itype(i))
  60 c       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
  61 c     &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
  62         phi(i+1)=gen_phi(i+1,it1,it)
  63         if (back) then
  64           phi(i)=gen_phi(i+1,it2,it1)
  65           print *,'phi(',i,')=',phi(i)
  66           theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
  67           if (it2.ne.10) then
  68             nsi=0
  69             fail=.true.
  70             do while (fail.and.nsi.le.maxsi)
  71               call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
  72               nsi=nsi+1
  73             enddo
  74             if (nsi.gt.maxsi) return1
  75           endif
  76           call locate_next_res(i-1)
  77         endif
  78         theta(i)=gen_theta(it1,phi(i),phi(i+1))
  79         if (it1.ne.10) then
  80         nsi=0
  81         fail=.true.
  82         do while (fail.and.nsi.le.maxsi)
  83           call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
  84           nsi=nsi+1
  85         enddo
  86         if (nsi.gt.maxsi) return1
  87         endif
  88         call locate_next_res(i)
  89         if (overlap(i-1)) then
  90           if (nit.lt.maxnit) then
  91             back=.true.
  92             nit=nit+1
  93           else
  94             nit=0
  95             if (i.gt.3) then
  96               back=.true.
  97               i=i-1
  98             else
  99               write (iout,'(a)')
 100      &  'Cannot generate non-overlaping conformation. Increase MAXNIT.'
 101               write (*,'(a)')
 102      &  'Cannot generate non-overlaping conformation. Increase MAXNIT.'
 103               return1
 104             endif
 105           endif
 106         else
 107           back=.false.
 108           nit=0
 109           i=i+1
 110         endif
 111         niter=niter+1
 112       enddo
 113       if (niter.ge.maxgen) then
 114         write (iout,'(a,2i5)')
 115      & 'Too many trials in conformation generation',niter,maxgen
 116         write (*,'(a,2i5)')
 117      & 'Too many trials in conformation generation',niter,maxgen
 118         return1
 119       endif
 120       do j=1,3
 121         c(j,nres+1)=c(j,1)
 122         c(j,nres+nres)=c(j,nres)
 123       enddo
 124       return
 125       end
 126 c-------------------------------------------------------------------------
 127       logical function overlap(i)
 128       implicit real*8 (a-h,o-z)
 129       include 'DIMENSIONS'
 130       include 'COMMON.CHAIN'
 131       include 'COMMON.INTERACT'
 132       include 'COMMON.FFIELD'
 133       data redfac /0.5D0/
 134       overlap=.false.
 135       iti=abs(itype(i))
 136       if (iti.gt.ntyp) return
 137 C Check for SC-SC overlaps.
 138 cd    print *,'nnt=',nnt,' nct=',nct
 139       do j=nnt,i-1
 140         itj=abs(itype(j))
 141         if (j.lt.i-1 .or. ipot.ne.4) then
 142           rcomp=sigmaii(iti,itj)
 143         else
 144           rcomp=sigma(iti,itj)
 145         endif
 146 cd      print *,'j=',j
 147         if (dist(nres+i,nres+j).lt.redfac*rcomp) then
 148           overlap=.true.
 149 c        print *,'overlap, SC-SC: i=',i,' j=',j,
 150 c     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
 151 c     &     rcomp
 152           return
 153         endif
 154       enddo
 155 C Check for overlaps between the added peptide group and the preceding
 156 C SCs.
 157       iteli=itel(i)
 158       do j=1,3
 159         c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
 160       enddo
 161       do j=nnt,i-2
 162         itj=abs(itype(j))
 163 cd      print *,'overlap, p-Sc: i=',i,' j=',j,
 164 cd   &         ' dist=',dist(nres+j,maxres2+1)
 165         if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
 166           overlap=.true.
 167           return
 168         endif
 169       enddo
 170 C Check for overlaps between the added side chain and the preceding peptide
 171 C groups.
 172       do j=1,nnt-2
 173         do k=1,3
 174           c(k,maxres2+1)=0.5D0*(c(k,j)+c(k,j+1))
 175         enddo
 176 cd      print *,'overlap, SC-p: i=',i,' j=',j,
 177 cd   &         ' dist=',dist(nres+i,maxres2+1)
 178         if (dist(nres+i,maxres2+1).lt.4.0D0*redfac) then
 179           overlap=.true.
 180           return
 181         endif
 182       enddo
 183 C Check for p-p overlaps
 184       do j=1,3
 185         c(j,maxres2+2)=0.5D0*(c(j,i)+c(j,i+1))
 186       enddo
 187       do j=nnt,i-2
 188         itelj=itel(j)
 189         do k=1,3
 190           c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
 191         enddo
 192 cd      print *,'overlap, p-p: i=',i,' j=',j,
 193 cd   &         ' dist=',dist(maxres2+1,maxres2+2)
 194         if(iteli.ne.0.and.itelj.ne.0)then
 195         if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfac) then
 196           overlap=.true.
 197           return
 198         endif
 199         endif
 200       enddo
 201       return
 202       end
 203 c--------------------------------------------------------------------------
 204       double precision function gen_phi(i,it1,it2)
 205       implicit real*8 (a-h,o-z)
 206       include 'DIMENSIONS'
 207       include 'COMMON.GEO'
 208       include 'COMMON.BOUNDS'
 209 c      gen_phi=ran_number(-pi,pi)
 210 C 8/13/98 Generate phi using pre-defined boundaries
 211       gen_phi=ran_number(phibound(1,i),phibound(2,i))
 212       return
 213       end
 214 c---------------------------------------------------------------------------
 215       double precision function gen_theta(it,gama,gama1)
 216       implicit real*8 (a-h,o-z)
 217       include 'DIMENSIONS'
 218       include 'COMMON.LOCAL'
 219       include 'COMMON.GEO'
 220       double precision y(2),z(2)
 221       double precision theta_max,theta_min
 222 c     print *,'gen_theta: it=',it
 223       theta_min=0.05D0*pi
 224       theta_max=0.95D0*pi
 225       if (dabs(gama).gt.dwapi) then
 226         y(1)=dcos(gama)
 227         y(2)=dsin(gama)
 228       else
 229         y(1)=0.0D0
 230         y(2)=0.0D0
 231       endif
 232       if (dabs(gama1).gt.dwapi) then
 233         z(1)=dcos(gama1)
 234         z(2)=dsin(gama1)
 235       else
 236         z(1)=0.0D0
 237         z(2)=0.0D0
 238       endif
 239       thet_pred_mean=a0thet(it)
 240       do k=1,2
 241         thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k)
 242      &                +bthet(k,it,1,1)*z(k)
 243       enddo
 244       sig=polthet(3,it)
 245       do j=2,0,-1
 246         sig=sig*thet_pred_mean+polthet(j,it)
 247       enddo
 248       sig=0.5D0/(sig*sig+sigc0(it))
 249       ak=dexp(gthet(1,it)-
 250      &0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
 251 c     print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
 252 c     print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
 253       theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak)
 254       if (theta_temp.lt.theta_min) theta_temp=theta_min
 255       if (theta_temp.gt.theta_max) theta_temp=theta_max
 256       gen_theta=theta_temp
 257 c     print '(a)','Exiting GENTHETA.'
 258       return
 259       end
 260 c-------------------------------------------------------------------------
 261       subroutine gen_side(it,the,al,om,fail)
 262       implicit real*8 (a-h,o-z)
 263       include 'DIMENSIONS'
 264       include 'COMMON.GEO'
 265       include 'COMMON.LOCAL'
 266       include 'COMMON.SETUP'
 267       include 'COMMON.IOUNITS'
 268       double precision MaxBoxLen /10.0D0/
 269       double precision Ap_inv(3,3),a(3,3),z(3,maxlob),W1(maxlob),
 270      & sumW(0:maxlob),y(2),cm(2),eig(2),box(2,2),work(100),detAp(maxlob)
 271       double precision eig_limit /1.0D-8/
 272       double precision Big /10.0D0/
 273       double precision vec(3,3)
 274       logical lprint,fail,lcheck
 275       lcheck=.false.
 276       lprint=.false.
 277       fail=.false.
 278       if (the.eq.0.0D0 .or. the.eq.pi) then
 279 #ifdef MPI
 280         write (*,'(a,i4,a,i3,a,1pe14.5)')
 281      & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
 282 #else
 283 cd        write (iout,'(a,i3,a,1pe14.5)')
 284 cd     &   'Error in GenSide: it=',it,' theta=',the
 285 #endif
 286         fail=.true.
 287         return
 288       endif
 289       tant=dtan(the-pipol)
 290       nlobit=nlob(it)
 291       if (lprint) then
 292 #ifdef MPI
 293         print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
 294         write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
 295 #endif
 296         print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
 297         write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,
 298      &     ' tant=',tant
 299       endif
 300       do i=1,nlobit
 301        zz1=tant-censc(1,i,it)
 302         do k=1,3
 303           do l=1,3
 304             a(k,l)=gaussc(k,l,i,it)
 305           enddo
 306         enddo
 307         detApi=a(2,2)*a(3,3)-a(2,3)**2
 308         Ap_inv(2,2)=a(3,3)/detApi
 309         Ap_inv(2,3)=-a(2,3)/detApi
 310         Ap_inv(3,2)=Ap_inv(2,3)
 311         Ap_inv(3,3)=a(2,2)/detApi
 312         if (lprint) then
 313           write (*,'(/a,i2/)') 'Cluster #',i
 314           write (*,'(3(1pe14.5),5x,1pe14.5)')
 315      &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
 316           write (iout,'(/a,i2/)') 'Cluster #',i
 317           write (iout,'(3(1pe14.5),5x,1pe14.5)')
 318      &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
 319         endif
 320         W1i=0.0D0
 321         do k=2,3
 322           do l=2,3
 323             W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
 324           enddo
 325         enddo
 326         W1i=a(1,1)-W1i
 327         W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
 328 c        if (lprint) write(*,'(a,3(1pe15.5)/)')
 329 c     &          'detAp, W1, anormi',detApi,W1i,anormi
 330         do k=2,3
 331           zk=censc(k,i,it)
 332           do l=2,3
 333             zk=zk+zz1*Ap_inv(k,l)*a(l,1)
 334           enddo
 335           z(k,i)=zk
 336         enddo
 337         detAp(i)=dsqrt(detApi)
 338       enddo
 339
 340       if (lprint) then
 341         print *,'W1:',(w1(i),i=1,nlobit)
 342         print *,'detAp:',(detAp(i),i=1,nlobit)
 343         print *,'Z'
 344         do i=1,nlobit
 345           print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
 346         enddo
 347         write (iout,*) 'W1:',(w1(i),i=1,nlobit)
 348         write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
 349         write (iout,*) 'Z'
 350         do i=1,nlobit
 351           write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
 352         enddo
 353       endif
 354       if (lcheck) then
 355 C Writing the distribution just to check the procedure
 356       fac=0.0D0
 357       dV=deg2rad**2*10.0D0
 358       sum=0.0D0
 359       sum1=0.0D0
 360       do i=1,nlobit
 361         fac=fac+W1(i)/detAp(i)
 362       enddo
 363       fac=1.0D0/(2.0D0*fac*pi)
 364 cd    print *,it,'fac=',fac
 365       do ial=90,180,2
 366         y(1)=deg2rad*ial
 367         do iom=-180,180,5
 368           y(2)=deg2rad*iom
 369           wart=0.0D0
 370           do i=1,nlobit
 371             do j=2,3
 372               do k=2,3
 373                 a(j-1,k-1)=gaussc(j,k,i,it)
 374               enddo
 375             enddo
 376             y2=y(2)
 377
 378             do iii=-1,1
 379
 380               y(2)=y2+iii*dwapi
 381
 382               wykl=0.0D0
 383               do j=1,2
 384                 do k=1,2
 385                   wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
 386                 enddo
 387               enddo
 388               wart=wart+W1(i)*dexp(-0.5D0*wykl)
 389
 390             enddo
 391
 392             y(2)=y2
 393
 394           enddo
 395 c         print *,'y',y(1),y(2),' fac=',fac
 396           wart=fac*wart
 397           write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
 398           sum=sum+wart
 399           sum1=sum1+1.0D0
 400         enddo
 401       enddo
 402 c     print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
 403       return
 404       endif
 405
 406 C Calculate the CM of the system
 407 C
 408       do i=1,nlobit
 409         W1(i)=W1(i)/detAp(i)
 410       enddo
 411       sumW(0)=0.0D0
 412       do i=1,nlobit
 413         sumW(i)=sumW(i-1)+W1(i)
 414       enddo
 415       cm(1)=z(2,1)*W1(1)
 416       cm(2)=z(3,1)*W1(1)
 417       do j=2,nlobit
 418         cm(1)=cm(1)+z(2,j)*W1(j)
 419         cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
 420       enddo
 421       cm(1)=cm(1)/sumW(nlobit)
 422       cm(2)=cm(2)/sumW(nlobit)
 423       if (cm(1).gt.Big .or. cm(1).lt.-Big .or.
 424      & cm(2).gt.Big .or. cm(2).lt.-Big) then
 425 cd        write (iout,'(a)')
 426 cd     & 'Unexpected error in GenSide - CM coordinates too large.'
 427 cd        write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
 428 cd        write (*,'(a)')
 429 cd     & 'Unexpected error in GenSide - CM coordinates too large.'
 430 cd        write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
 431         fail=.true.
 432         return
 433       endif
 434 cd    print *,'CM:',cm(1),cm(2)
 435 C
 436 C Find the largest search distance from CM
 437 C
 438       radmax=0.0D0
 439       do i=1,nlobit
 440         do j=2,3
 441           do k=2,3
 442             a(j-1,k-1)=gaussc(j,k,i,it)
 443           enddo
 444         enddo
 445 #ifdef NAG
 446         call f02faf('N','U',2,a,3,eig,work,100,ifail)
 447 #else
 448         call djacob(2,3,10000,1.0d-10,a,vec,eig)
 449 #endif
 450 #ifdef MPI
 451         if (lprint) then
 452           print *,'*************** CG Processor',me
 453           print *,'CM:',cm(1),cm(2)
 454           write (iout,*) '*************** CG Processor',me
 455           write (iout,*) 'CM:',cm(1),cm(2)
 456           print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
 457           write (iout,'(A,8f10.5)')
 458      &        'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
 459         endif
 460 #endif
 461         if (eig(1).lt.eig_limit) then
 462           write(iout,'(a)')
 463      &     'From Mult_Norm: Eigenvalues of A are too small.'
 464           write(*,'(a)')
 465      &     'From Mult_Norm: Eigenvalues of A are too small.'
 466           fail=.true.
 467           return
 468         endif
 469         radius=0.0D0
 470 cd      print *,'i=',i
 471         do j=1,2
 472           radius=radius+pinorm(z(j+1,i)-cm(j))**2
 473         enddo
 474         radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
 475         if (radius.gt.radmax) radmax=radius
 476       enddo
 477       if (radmax.gt.pi) radmax=pi
 478 C
 479 C Determine the boundaries of the search rectangle.
 480 C
 481       if (lprint) then
 482         print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
 483         print '(a,4(1pe14.4))','radmax: ',radmax
 484       endif
 485       box(1,1)=dmax1(cm(1)-radmax,0.0D0)
 486       box(2,1)=dmin1(cm(1)+radmax,pi)
 487       box(1,2)=cm(2)-radmax
 488       box(2,2)=cm(2)+radmax
 489       if (lprint) then
 490 #ifdef MPI
 491         print *,'CG Processor',me,' Array BOX:'
 492 #else
 493         print *,'Array BOX:'
 494 #endif
 495         print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
 496         print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
 497 #ifdef MPI
 498         write (iout,*)'CG Processor',me,' Array BOX:'
 499 #else
 500         write (iout,*)'Array BOX:'
 501 #endif
 502         write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
 503         write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
 504       endif
 505       if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
 506 #ifdef MPI
 507         write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
 508         write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
 509 #else
 510 c        write (iout,'(a)') 'Bad sampling box.'
 511 #endif
 512         fail=.true.
 513         return
 514       endif
 515       which_lobe=ran_number(0.0D0,sumW(nlobit))
 516 c     print '(a,1pe14.4)','which_lobe=',which_lobe
 517       do i=1,nlobit
 518         if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
 519       enddo
 520     1 ilob=i
 521 c     print *,'ilob=',ilob,' nlob=',nlob(it)
 522       do i=2,3
 523         cm(i-1)=z(i,ilob)
 524         do j=2,3
 525           a(i-1,j-1)=gaussc(i,j,ilob,it)
 526         enddo
 527       enddo
 528 cd    print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
 529       call mult_norm1(3,2,a,cm,box,y,fail)
 530       if (fail) return
 531       al=y(1)
 532       om=pinorm(y(2))
 533 cd    print *,'al=',al,' om=',om
 534 cd    stop
 535       return
 536       end
 537 c---------------------------------------------------------------------------
 538       double precision function ran_number(x1,x2)
 539 C Calculate a random real number from the range (x1,x2).
 540       implicit real*8 (a-h,o-z)
 541       include 'DIMENSIONS'
 542       double precision x1,x2,fctor
 543       data fctor /2147483647.0D0/
 544 #ifdef MPI
 545       include "mpif.h"
 546       include 'COMMON.SETUP'
 547       ran_number=x1+(x2-x1)*prng_next(me)
 548 #else
 549       call vrnd(ix,1)
 550       ran_number=x1+(x2-x1)*ix/fctor
 551 #endif
 552       return
 553       end
 554 c--------------------------------------------------------------------------
 555       integer function iran_num(n1,n2)
 556 C Calculate a random integer number from the range (n1,n2).
 557       implicit real*8 (a-h,o-z)
 558       include 'DIMENSIONS'
 559       integer n1,n2,ix
 560       real fctor /2147483647.0/
 561 #ifdef MPI
 562       include "mpif.h"
 563       include 'COMMON.SETUP'
 564       ix=n1+(n2-n1+1)*prng_next(me)
 565       if (ix.lt.n1) ix=n1
 566       if (ix.gt.n2) ix=n2
 567       iran_num=ix
 568 #else
 569       call vrnd(ix,1)
 570       ix=n1+(n2-n1+1)*(ix/fctor)
 571       if (ix.gt.n2) ix=n2
 572       iran_num=ix
 573 #endif
 574       return
 575       end
 576 c--------------------------------------------------------------------------
 577       double precision function binorm(x1,x2,sigma1,sigma2,ak)
 578       implicit real*8 (a-h,o-z)
 579 c     print '(a)','Enter BINORM.'
 580       alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
 581       aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
 582       seg=sigma1/(sigma1+ak*sigma2)
 583       alen=ran_number(0.0D0,1.0D0)
 584       if (alen.lt.seg) then
 585         binorm=anorm_distr(x1,sigma1,alowb,aupb)
 586       else
 587         binorm=anorm_distr(x2,sigma2,alowb,aupb)
 588       endif
 589 c     print '(a)','Exiting BINORM.'
 590       return
 591       end
 592 c-----------------------------------------------------------------------
 593 c      double precision function anorm_distr(x,sigma,alowb,aupb)
 594 c      implicit real*8 (a-h,o-z)
 595 c     print '(a)','Enter ANORM_DISTR.'
 596 c   10 y=ran_number(alowb,aupb)
 597 c      expon=dexp(-0.5D0*((y-x)/sigma)**2)
 598 c      ran=ran_number(0.0D0,1.0D0)
 599 c      if (expon.lt.ran) goto 10
 600 c      anorm_distr=y
 601 c     print '(a)','Exiting ANORM_DISTR.'
 602 c      return
 603 c      end
 604 c-----------------------------------------------------------------------
 605         double precision function anorm_distr(x,sigma,alowb,aupb)
 606         implicit real*8 (a-h,o-z)
 607 c  to make a normally distributed deviate with zero mean and unit variance
 608 c
 609         integer iset
 610         real fac,gset,rsq,v1,v2,ran1
 611         save iset,gset
 612         data iset/0/
 613         if(iset.eq.0) then
 614 1               v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
 615                 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
 616                 rsq=v1**2+v2**2
 617                 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
 618                 fac=sqrt(-2.0d0*log(rsq)/rsq)
 619                 gset=v1*fac
 620                 gaussdev=v2*fac
 621                 iset=1
 622         else
 623                 gaussdev=gset
 624                 iset=0
 625         endif
 626         anorm_distr=x+gaussdev*sigma
 627       return
 628       end
 629 c------------------------------------------------------------------------
 630       subroutine mult_norm(lda,n,a,x,fail)
 631 C
 632 C Generate the vector X whose elements obey the multiple-normal distribution
 633 C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
 634 C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
 635 C eigenvalue of the matrix A is close to 0.
 636 C
 637       implicit double precision (a-h,o-z)
 638       double precision a(lda,n),x(n),eig(100),vec(3,3),work(100)
 639       double precision eig_limit /1.0D-8/
 640       logical fail
 641       fail=.false.
 642 c     print '(a)','Enter MULT_NORM.'
 643 C
 644 C Find the smallest eigenvalue of the matrix A.
 645 C
 646 c     do i=1,n
 647 c       print '(8f10.5)',(a(i,j),j=1,n)
 648 c     enddo
 649 #ifdef NAG
 650       call f02faf('V','U',2,a,lda,eig,work,100,ifail)
 651 #else
 652       call djacob(2,lda,10000,1.0d-10,a,vec,eig)
 653 #endif
 654 c     print '(8f10.5)',(eig(i),i=1,n)
 655 C     print '(a)'
 656 c     do i=1,n
 657 c       print '(8f10.5)',(a(i,j),j=1,n)
 658 c     enddo
 659       if (eig(1).lt.eig_limit) then
 660         print *,'From Mult_Norm: Eigenvalues of A are too small.'
 661         fail=.true.
 662         return
 663       endif
 664 C
 665 C Generate points following the normal distributions along the principal
 666 C axes of the moment matrix. Store in WORK.
 667 C
 668       do i=1,n
 669         sigma=1.0D0/dsqrt(eig(i))
 670         alim=-3.0D0*sigma
 671         work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
 672       enddo
 673 C
 674 C Transform the vector of normal variables back to the original basis.
 675 C
 676       do i=1,n
 677         xi=0.0D0
 678         do j=1,n
 679           xi=xi+a(i,j)*work(j)
 680         enddo
 681         x(i)=xi
 682       enddo
 683       return
 684       end
 685 c------------------------------------------------------------------------
 686       subroutine mult_norm1(lda,n,a,z,box,x,fail)
 687 C
 688 C Generate the vector X whose elements obey the multi-gaussian multi-dimensional
 689 C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
 690 C leading dimension of the moment matrix A, n is the dimension of the
 691 C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
 692 C smallest eigenvalue of the matrix A is close to 0.
 693 C
 694       implicit real*8 (a-h,o-z)
 695       include 'DIMENSIONS'
 696 #ifdef MPI
 697       include 'mpif.h'
 698 #endif
 699       double precision a(lda,n),z(n),x(n),box(n,n)
 700       double precision etmp
 701       include 'COMMON.IOUNITS'
 702 #ifdef MP
 703       include 'COMMON.SETUP'
 704 #endif
 705       logical fail
 706 C
 707 C Generate points following the normal distributions along the principal
 708 C axes of the moment matrix. Store in WORK.
 709 C
 710 cd    print *,'CG Processor',me,' entered MultNorm1.'
 711 cd    print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
 712 cd    do i=1,n
 713 cd      print *,i,box(1,i),box(2,i)
 714 cd    enddo
 715       istep = 0
 716    10 istep = istep + 1
 717       if (istep.gt.10000) then
 718 c        write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
 719 c     & ' in MultNorm1.'
 720 c        write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
 721 c     & ' in MultNorm1.'
 722 c        write (iout,*) 'box',box
 723 c        write (iout,*) 'a',a
 724 c        write (iout,*) 'z',z
 725         fail=.true.
 726         return
 727       endif
 728       do i=1,n
 729        x(i)=ran_number(box(1,i),box(2,i))
 730       enddo
 731       ww=0.0D0
 732       do i=1,n
 733         xi=pinorm(x(i)-z(i))
 734         ww=ww+0.5D0*a(i,i)*xi*xi
 735         do j=i+1,n
 736           ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
 737         enddo
 738       enddo
 739       dec=ran_number(0.0D0,1.0D0)
 740 c      print *,(x(i),i=1,n),ww,dexp(-ww),dec
 741 crc   if (dec.gt.dexp(-ww)) goto 10
 742       if(-ww.lt.100) then
 743        etmp=dexp(-ww)
 744       else
 745        return
 746       endif
 747       if (dec.gt.etmp) goto 10
 748 cd    print *,'CG Processor',me,' exitting MultNorm1.'
 749       return
 750       end
 751 c
 752 crc--------------------------------------
 753       subroutine overlap_sc(scfail)
 754 c     Internal and cartesian coordinates must be consistent as input,
 755 c     and will be up-to-date on return.
 756 c     At the end of this procedure, scfail is true if there are
 757 c     overlapping residues left, or false otherwise (success)
 758       implicit real*8 (a-h,o-z)
 759       include 'DIMENSIONS'
 760       include 'COMMON.CHAIN'
 761       include 'COMMON.INTERACT'
 762       include 'COMMON.FFIELD'
 763       include 'COMMON.VAR'
 764       include 'COMMON.SBRIDGE'
 765       include 'COMMON.IOUNITS'
 766       logical had_overlaps,fail,scfail
 767       integer ioverlap(maxres),ioverlap_last
 768
 769       had_overlaps=.false.
 770       call overlap_sc_list(ioverlap,ioverlap_last)
 771       if (ioverlap_last.gt.0) then
 772         write (iout,*) '#OVERLAPing residues ',ioverlap_last
 773         write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
 774         had_overlaps=.true.
 775       endif
 776
 777       maxsi=1000
 778       do k=1,1000
 779         if (ioverlap_last.eq.0) exit
 780
 781         do ires=1,ioverlap_last
 782           i=ioverlap(ires)
 783           iti=abs(itype(i))
 784           if (iti.ne.10) then
 785             nsi=0
 786             fail=.true.
 787             do while (fail.and.nsi.le.maxsi)
 788               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 789               nsi=nsi+1
 790             enddo
 791             if(fail) goto 999
 792           endif
 793         enddo
 794
 795         call chainbuild
 796         call overlap_sc_list(ioverlap,ioverlap_last)
 797 c        write (iout,*) 'Overlaping residues ',ioverlap_last,
 798 c     &           (ioverlap(j),j=1,ioverlap_last)
 799       enddo
 800
 801       if (k.le.1000.and.ioverlap_last.eq.0) then
 802         scfail=.false.
 803         if (had_overlaps) then
 804           write (iout,*) '#OVERLAPing all corrected after ',k,
 805      &         ' random generation'
 806         endif
 807       else
 808         scfail=.true.
 809         write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
 810         write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
 811       endif
 812
 813       return
 814
 815  999  continue
 816       write (iout,'(a30,i5,a12,i4)')
 817      &               '#OVERLAP FAIL in gen_side after',maxsi,
 818      &               'iter for RES',i
 819       scfail=.true.
 820       return
 821       end
 822
 823       subroutine overlap_sc_list(ioverlap,ioverlap_last)
 824       implicit real*8 (a-h,o-z)
 825       include 'DIMENSIONS'
 826       include 'COMMON.GEO'
 827       include 'COMMON.LOCAL'
 828       include 'COMMON.IOUNITS'
 829       include 'COMMON.CHAIN'
 830       include 'COMMON.INTERACT'
 831       include 'COMMON.FFIELD'
 832       include 'COMMON.VAR'
 833       include 'COMMON.CALC'
 834       logical fail
 835       integer ioverlap(maxres),ioverlap_last
 836       data redfac /0.5D0/
 837
 838       ioverlap_last=0
 839 C Check for SC-SC overlaps and mark residues
 840 c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
 841       ind=0
 842       do i=iatsc_s,iatsc_e
 843         itypi=abs(itype(i))
 844         itypi1=abs(itype(i+1))
 845         xi=c(1,nres+i)
 846         yi=c(2,nres+i)
 847         zi=c(3,nres+i)
 848         dxi=dc_norm(1,nres+i)
 849         dyi=dc_norm(2,nres+i)
 850         dzi=dc_norm(3,nres+i)
 851         dsci_inv=dsc_inv(itypi)
 852 c
 853        do iint=1,nint_gr(i)
 854          do j=istart(i,iint),iend(i,iint)
 855             ind=ind+1
 856             itypj=itype(j)
 857             dscj_inv=dsc_inv(itypj)
 858             sig0ij=sigma(itypi,itypj)
 859             chi1=chi(itypi,itypj)
 860             chi2=chi(itypj,itypi)
 861             chi12=chi1*chi2
 862             chip1=chip(itypi)
 863             chip2=chip(itypj)
 864             chip12=chip1*chip2
 865             alf1=alp(itypi)
 866             alf2=alp(itypj)
 867             alf12=0.5D0*(alf1+alf2)
 868           if (j.gt.i+1) then
 869            rcomp=sigmaii(itypi,itypj)
 870           else
 871            rcomp=sigma(itypi,itypj)
 872           endif
 873 c         print '(2(a3,2i3),a3,2f10.5)',
 874 c     &        ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
 875 c     &        ,rcomp
 876             xj=c(1,nres+j)-xi
 877             yj=c(2,nres+j)-yi
 878             zj=c(3,nres+j)-zi
 879             dxj=dc_norm(1,nres+j)
 880             dyj=dc_norm(2,nres+j)
 881             dzj=dc_norm(3,nres+j)
 882             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 883             rij=dsqrt(rrij)
 884             call sc_angular
 885             sigsq=1.0D0/sigsq
 886             sig=sig0ij*dsqrt(sigsq)
 887             rij_shift=1.0D0/rij-sig+sig0ij
 888
 889 ct          if ( 1.0/rij .lt. redfac*rcomp .or.
 890 ct     &       rij_shift.le.0.0D0 ) then
 891             if ( rij_shift.le.0.0D0 ) then
 892 cd           write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
 893 cd     &     'overlap SC-SC: i=',i,' j=',j,
 894 cd     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
 895 cd     &     rcomp,1.0/rij,rij_shift
 896           ioverlap_last=ioverlap_last+1
 897           ioverlap(ioverlap_last)=i
 898           do k=1,ioverlap_last-1
 899            if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
 900           enddo
 901           ioverlap_last=ioverlap_last+1
 902           ioverlap(ioverlap_last)=j
 903           do k=1,ioverlap_last-1
 904            if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
 905           enddo
 906          endif
 907         enddo
 908        enddo
 909       enddo
 910       return
 911       end