source/unres/src-HCD-5D/gen_rand_conf.F

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