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,lprn /.false./
 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         if (lprn) then
 518         write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
 519         write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
 520         endif
 521 #else
 522         if (lprn) write (iout,'(a)') 'Bad sampling box.'
 523 #endif
 524         fail=.true.
 525         return
 526       endif
 527       which_lobe=ran_number(0.0D0,sumW(nlobit))
 528 c     print '(a,1pe14.4)','which_lobe=',which_lobe
 529       do i=1,nlobit
 530         if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
 531       enddo
 532     1 ilob=i
 533 c     print *,'ilob=',ilob,' nlob=',nlob(it)
 534       do i=2,3
 535         cm(i-1)=z(i,ilob)
 536         do j=2,3
 537           a(i-1,j-1)=gaussc(i,j,ilob,it)
 538         enddo
 539       enddo
 540 cd    print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
 541       call mult_norm1(3,2,a,cm,box,y,fail)
 542       if (fail) return
 543       al=y(1)
 544       om=pinorm(y(2))
 545 cd    print *,'al=',al,' om=',om
 546 cd    stop
 547       return
 548       end
 549 c---------------------------------------------------------------------------
 550       double precision function ran_number(x1,x2)
 551 C Calculate a random real number from the range (x1,x2).
 552       implicit real*8 (a-h,o-z)
 553       include 'DIMENSIONS'
 554       double precision x1,x2,fctor
 555       data fctor /2147483647.0D0/
 556 #ifdef MPI
 557       include "mpif.h"
 558       include 'COMMON.SETUP'
 559       ran_number=x1+(x2-x1)*prng_next(me)
 560 #else
 561       call vrnd(ix,1)
 562       ran_number=x1+(x2-x1)*ix/fctor
 563 #endif
 564       return
 565       end
 566 c--------------------------------------------------------------------------
 567       integer function iran_num(n1,n2)
 568 C Calculate a random integer number from the range (n1,n2).
 569       implicit real*8 (a-h,o-z)
 570       include 'DIMENSIONS'
 571       integer n1,n2,ix
 572       real fctor /2147483647.0/
 573 #ifdef MPI
 574       include "mpif.h"
 575       include 'COMMON.SETUP'
 576       ix=n1+(n2-n1+1)*prng_next(me)
 577       if (ix.lt.n1) ix=n1
 578       if (ix.gt.n2) ix=n2
 579       iran_num=ix
 580 #else
 581       call vrnd(ix,1)
 582       ix=n1+(n2-n1+1)*(ix/fctor)
 583       if (ix.gt.n2) ix=n2
 584       iran_num=ix
 585 #endif
 586       return
 587       end
 588 c--------------------------------------------------------------------------
 589       double precision function binorm(x1,x2,sigma1,sigma2,ak)
 590       implicit real*8 (a-h,o-z)
 591 c     print '(a)','Enter BINORM.'
 592       alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
 593       aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
 594       seg=sigma1/(sigma1+ak*sigma2)
 595       alen=ran_number(0.0D0,1.0D0)
 596       if (alen.lt.seg) then
 597         binorm=anorm_distr(x1,sigma1,alowb,aupb)
 598       else
 599         binorm=anorm_distr(x2,sigma2,alowb,aupb)
 600       endif
 601 c     print '(a)','Exiting BINORM.'
 602       return
 603       end
 604 c-----------------------------------------------------------------------
 605 c      double precision function anorm_distr(x,sigma,alowb,aupb)
 606 c      implicit real*8 (a-h,o-z)
 607 c     print '(a)','Enter ANORM_DISTR.'
 608 c   10 y=ran_number(alowb,aupb)
 609 c      expon=dexp(-0.5D0*((y-x)/sigma)**2)
 610 c      ran=ran_number(0.0D0,1.0D0)
 611 c      if (expon.lt.ran) goto 10
 612 c      anorm_distr=y
 613 c     print '(a)','Exiting ANORM_DISTR.'
 614 c      return
 615 c      end
 616 c-----------------------------------------------------------------------
 617         double precision function anorm_distr(x,sigma,alowb,aupb)
 618         implicit real*8 (a-h,o-z)
 619 c  to make a normally distributed deviate with zero mean and unit variance
 620 c
 621         integer iset
 622         real fac,gset,rsq,v1,v2,ran1
 623         save iset,gset
 624         data iset/0/
 625         if(iset.eq.0) then
 626 1               v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
 627                 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
 628                 rsq=v1**2+v2**2
 629                 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
 630                 fac=sqrt(-2.0d0*log(rsq)/rsq)
 631                 gset=v1*fac
 632                 gaussdev=v2*fac
 633                 iset=1
 634         else
 635                 gaussdev=gset
 636                 iset=0
 637         endif
 638         anorm_distr=x+gaussdev*sigma
 639         return
 640         end
 641 c------------------------------------------------------------------------
 642       subroutine mult_norm(lda,n,a,x,fail)
 643 C
 644 C Generate the vector X whose elements obey the multiple-normal distribution
 645 C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
 646 C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
 647 C eigenvalue of the matrix A is close to 0.
 648 C
 649       implicit double precision (a-h,o-z)
 650       double precision a(lda,n),x(n),eig(100),vec(3,3),work(100)
 651       double precision eig_limit /1.0D-8/
 652       logical fail
 653       fail=.false.
 654 c     print '(a)','Enter MULT_NORM.'
 655 C
 656 C Find the smallest eigenvalue of the matrix A.
 657 C
 658 c     do i=1,n
 659 c       print '(8f10.5)',(a(i,j),j=1,n)
 660 c     enddo
 661 #ifdef NAG
 662       call f02faf('V','U',2,a,lda,eig,work,100,ifail)
 663 #else
 664       call djacob(2,lda,10000,1.0d-10,a,vec,eig)
 665 #endif
 666 c     print '(8f10.5)',(eig(i),i=1,n)
 667 C     print '(a)'
 668 c     do i=1,n
 669 c       print '(8f10.5)',(a(i,j),j=1,n)
 670 c     enddo
 671       if (eig(1).lt.eig_limit) then
 672         print *,'From Mult_Norm: Eigenvalues of A are too small.'
 673         fail=.true.
 674         return
 675       endif
 676 C
 677 C Generate points following the normal distributions along the principal
 678 C axes of the moment matrix. Store in WORK.
 679 C
 680       do i=1,n
 681         sigma=1.0D0/dsqrt(eig(i))
 682         alim=-3.0D0*sigma
 683         work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
 684       enddo
 685 C
 686 C Transform the vector of normal variables back to the original basis.
 687 C
 688       do i=1,n
 689         xi=0.0D0
 690         do j=1,n
 691           xi=xi+a(i,j)*work(j)
 692         enddo
 693         x(i)=xi
 694       enddo
 695       return
 696       end
 697 c------------------------------------------------------------------------
 698       subroutine mult_norm1(lda,n,a,z,box,x,fail)
 699 C
 700 C Generate the vector X whose elements obey the multi-gaussian multi-dimensional
 701 C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
 702 C leading dimension of the moment matrix A, n is the dimension of the
 703 C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
 704 C smallest eigenvalue of the matrix A is close to 0.
 705 C
 706       implicit real*8 (a-h,o-z)
 707       include 'DIMENSIONS'
 708 #ifdef MPI
 709       include 'mpif.h'
 710 #endif
 711       double precision a(lda,n),z(n),x(n),box(n,n)
 712       double precision etmp
 713       include 'COMMON.IOUNITS'
 714 #ifdef MP
 715       include 'COMMON.SETUP'
 716 #endif
 717       logical fail
 718 C
 719 C Generate points following the normal distributions along the principal
 720 C axes of the moment matrix. Store in WORK.
 721 C
 722 cd    print *,'CG Processor',me,' entered MultNorm1.'
 723 cd    print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
 724 cd    do i=1,n
 725 cd      print *,i,box(1,i),box(2,i)
 726 cd    enddo
 727       istep = 0
 728    10 istep = istep + 1
 729       if (istep.gt.10000) then
 730 c        write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
 731 c     & ' in MultNorm1.'
 732 c        write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
 733 c     & ' in MultNorm1.'
 734 c        write (iout,*) 'box',box
 735 c        write (iout,*) 'a',a
 736 c        write (iout,*) 'z',z
 737         fail=.true.
 738         return
 739       endif
 740       do i=1,n
 741        x(i)=ran_number(box(1,i),box(2,i))
 742       enddo
 743       ww=0.0D0
 744       do i=1,n
 745         xi=pinorm(x(i)-z(i))
 746         ww=ww+0.5D0*a(i,i)*xi*xi
 747         do j=i+1,n
 748           ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
 749         enddo
 750       enddo
 751       dec=ran_number(0.0D0,1.0D0)
 752 c      print *,(x(i),i=1,n),ww,dexp(-ww),dec
 753 crc   if (dec.gt.dexp(-ww)) goto 10
 754       if(-ww.lt.100) then
 755        etmp=dexp(-ww)
 756       else
 757        return
 758       endif
 759       if (dec.gt.etmp) goto 10
 760 cd    print *,'CG Processor',me,' exitting MultNorm1.'
 761       return
 762       end
 763 c
 764 crc--------------------------------------
 765       subroutine overlap_sc(scfail)
 766 c     Internal and cartesian coordinates must be consistent as input,
 767 c     and will be up-to-date on return.
 768 c     At the end of this procedure, scfail is true if there are
 769 c     overlapping residues left, or false otherwise (success)
 770       implicit real*8 (a-h,o-z)
 771       include 'DIMENSIONS'
 772       include 'COMMON.CHAIN'
 773       include 'COMMON.INTERACT'
 774       include 'COMMON.FFIELD'
 775       include 'COMMON.VAR'
 776       include 'COMMON.SBRIDGE'
 777       include 'COMMON.IOUNITS'
 778       logical had_overlaps,fail,scfail
 779       integer ioverlap(maxres),ioverlap_last
 780
 781       had_overlaps=.false.
 782       call overlap_sc_list(ioverlap,ioverlap_last)
 783       if (ioverlap_last.gt.0) then
 784         write (iout,*) '#OVERLAPing residues ',ioverlap_last
 785         write (iout,'(18i5)') (ioverlap(k),k=1,ioverlap_last)
 786         had_overlaps=.true.
 787       endif
 788
 789       maxsi=1000
 790       do k=1,1000
 791         if (ioverlap_last.eq.0) exit
 792
 793         do ires=1,ioverlap_last
 794           i=ioverlap(ires)
 795           iti=iabs(itype(i))
 796           if (iti.ne.10) then
 797             nsi=0
 798             fail=.true.
 799             do while (fail.and.nsi.le.maxsi)
 800               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 801               nsi=nsi+1
 802             enddo
 803             if(fail) goto 999
 804           endif
 805         enddo
 806
 807         call chainbuild_extconf
 808         call overlap_sc_list(ioverlap,ioverlap_last)
 809         write (iout,*) 'Overlaping residues ',ioverlap_last,
 810      &           (ioverlap(j),j=1,ioverlap_last)
 811       enddo
 812
 813       if (k.le.1000.and.ioverlap_last.eq.0) then
 814         scfail=.false.
 815         if (had_overlaps) then
 816           write (iout,*) '#OVERLAPing all corrected after ',k,
 817      &         ' random generation'
 818         endif
 819       else
 820         scfail=.true.
 821         write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
 822         write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
 823       endif
 824
 825       return
 826
 827  999  continue
 828       write (iout,'(a30,i5,a12,i4)')
 829      &               '#OVERLAP FAIL in gen_side after',maxsi,
 830      &               'iter for RES',i
 831       scfail=.true.
 832       return
 833       end
 834
 835       subroutine overlap_sc_list(ioverlap,ioverlap_last)
 836       implicit real*8 (a-h,o-z)
 837       include 'DIMENSIONS'
 838       include 'COMMON.GEO'
 839       include 'COMMON.LOCAL'
 840       include 'COMMON.IOUNITS'
 841       include 'COMMON.CHAIN'
 842       include 'COMMON.INTERACT'
 843       include 'COMMON.FFIELD'
 844       include 'COMMON.VAR'
 845       include 'COMMON.CALC'
 846       logical fail
 847       integer ioverlap(maxres),ioverlap_last
 848       data redfac /0.5D0/
 849
 850       write (iout,*) "overlap_sc_list"
 851 c      write(iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
 852       write(iout,*) "nnt",nnt," nct",nct
 853       ioverlap_last=0
 854 C Check for SC-SC overlaps and mark residues
 855 c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
 856       ind=0
 857 c      do i=iatsc_s,iatsc_e
 858       do i=nnt,nct
 859         itypi=iabs(itype(i))
 860         itypi1=iabs(itype(i+1))
 861         if (itypi.eq.ntyp1) cycle
 862         xi=c(1,nres+i)
 863         yi=c(2,nres+i)
 864         zi=c(3,nres+i)
 865         dxi=dc_norm(1,nres+i)
 866         dyi=dc_norm(2,nres+i)
 867         dzi=dc_norm(3,nres+i)
 868         dsci_inv=dsc_inv(itypi)
 869 c
 870 c        do iint=1,nint_gr(i)
 871 c          do j=istart(i,iint),iend(i,iint)
 872           do j=i+1,nct
 873             ind=ind+1
 874             itypj=iabs(itype(j))
 875             if (itypj.eq.ntyp1) cycle
 876 c           write (iout,*) "i,j",i,j," itypi,itypj",itypi,itypj
 877             dscj_inv=dsc_inv(itypj)
 878             sig0ij=sigma(itypi,itypj)
 879             chi1=chi(itypi,itypj)
 880             chi2=chi(itypj,itypi)
 881             chi12=chi1*chi2
 882             chip1=chip(itypi)
 883             chip2=chip(itypj)
 884             chip12=chip1*chip2
 885             alf1=alp(itypi)
 886             alf2=alp(itypj)
 887             alf12=0.5D0*(alf1+alf2)
 888             if (j.gt.i+1) then
 889               rcomp=sigmaii(itypi,itypj)
 890             else
 891               rcomp=sigma(itypi,itypj)
 892             endif
 893 c            write (iout,'(2(a3,2i5),a3,2f10.5)'),
 894 c     &        ' i=',i,itypi,' j=',j,itypj,' d=',dist(nres+i,nres+j)
 895 c     &        ,rcomp
 896             xj=c(1,nres+j)-xi
 897             yj=c(2,nres+j)-yi
 898             zj=c(3,nres+j)-zi
 899             dxj=dc_norm(1,nres+j)
 900             dyj=dc_norm(2,nres+j)
 901             dzj=dc_norm(3,nres+j)
 902             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 903             rij=dsqrt(rrij)
 904             call sc_angular
 905 c            write (iout,*) "dxj",dxj," dyj",dyj," dzj",dzj
 906 c            write (iout,*) "erij",erij
 907 c            write (iout,*) "om1",om1," om2",om2," om12",om12,
 908 c     &       " faceps1",faceps1," eps1",eps1
 909 c            write (iout,*) "sigsq",sigsq
 910             sigsq=1.0D0/sigsq
 911             sig=sig0ij*dsqrt(sigsq)
 912             rij_shift=1.0D0/rij-sig+sig0ij
 913 c            write (iout,*) "rij",1.0d0/rij," rij_shift",rij_shift,
 914 c     &       " sig",sig," sig0ij",sig0ij
 915 c            if ( rij_shift.le.0.0D0 ) then
 916             if ( rij_shift/sig0ij.le.0.1D0 ) then
 917 c              write (iout,*) "overlap",i,j
 918               write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)')
 919      &         'overlap SC-SC: i=',i,' j=',j,
 920      &         ' dist=',dist(nres+i,nres+j),' rcomp=',
 921      &         rcomp,1.0/rij,rij_shift
 922               ioverlap_last=ioverlap_last+1
 923               ioverlap(ioverlap_last)=i
 924               do k=1,ioverlap_last-1
 925                 if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
 926               enddo
 927               ioverlap_last=ioverlap_last+1
 928               ioverlap(ioverlap_last)=j
 929               do k=1,ioverlap_last-1
 930                 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
 931               enddo
 932             endif
 933           enddo
 934 c        enddo
 935       enddo
 936       return
 937       end