source/unres/src_MIN/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=itype(2)
  19         phi(4)=gen_phi(4,itype(2),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=itype(i-1)
  58         it2=itype(i-2)
  59         it=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 c         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=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=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=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)*y(k)+bthet(k,it)*z(k)
 242       enddo
 243       sig=polthet(3,it)
 244       do j=2,0,-1
 245         sig=sig*thet_pred_mean+polthet(j,it)
 246       enddo
 247       sig=0.5D0/(sig*sig+sigc0(it))
 248       ak=dexp(gthet(1,it)-
 249      &0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
 250 c     print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
 251 c     print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
 252       theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak)
 253       if (theta_temp.lt.theta_min) theta_temp=theta_min
 254       if (theta_temp.gt.theta_max) theta_temp=theta_max
 255       gen_theta=theta_temp
 256 c     print '(a)','Exiting GENTHETA.'
 257       return
 258       end
 259 c-------------------------------------------------------------------------
 260       subroutine gen_side(it,the,al,om,fail)
 261       implicit real*8 (a-h,o-z)
 262       include 'DIMENSIONS'
 263       include 'COMMON.GEO'
 264       include 'COMMON.LOCAL'
 265       include 'COMMON.SETUP'
 266       include 'COMMON.IOUNITS'
 267       double precision MaxBoxLen /10.0D0/
 268       double precision Ap_inv(3,3),a(3,3),z(3,maxlob),W1(maxlob),
 269      & sumW(0:maxlob),y(2),cm(2),eig(2),box(2,2),work(100),detAp(maxlob)
 270       double precision eig_limit /1.0D-8/
 271       double precision Big /10.0D0/
 272       double precision vec(3,3)
 273       logical lprint,fail,lcheck
 274       lcheck=.false.
 275       lprint=.false.
 276       fail=.false.
 277       if (the.eq.0.0D0 .or. the.eq.pi) then
 278 #ifdef MPI
 279         write (*,'(a,i4,a,i3,a,1pe14.5)')
 280      & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
 281 #else
 282 cd        write (iout,'(a,i3,a,1pe14.5)')
 283 cd     &   'Error in GenSide: it=',it,' theta=',the
 284 #endif
 285         fail=.true.
 286         return
 287       endif
 288       tant=dtan(the-pipol)
 289       nlobit=nlob(it)
 290       if (lprint) then
 291 #ifdef MPI
 292         print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
 293         write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
 294 #endif
 295         print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
 296         write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,
 297      &     ' tant=',tant
 298       endif
 299       do i=1,nlobit
 300        zz1=tant-censc(1,i,it)
 301         do k=1,3
 302           do l=1,3
 303             a(k,l)=gaussc(k,l,i,it)
 304           enddo
 305         enddo
 306         detApi=a(2,2)*a(3,3)-a(2,3)**2
 307         Ap_inv(2,2)=a(3,3)/detApi
 308         Ap_inv(2,3)=-a(2,3)/detApi
 309         Ap_inv(3,2)=Ap_inv(2,3)
 310         Ap_inv(3,3)=a(2,2)/detApi
 311         if (lprint) then
 312           write (*,'(/a,i2/)') 'Cluster #',i
 313           write (*,'(3(1pe14.5),5x,1pe14.5)')
 314      &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
 315           write (iout,'(/a,i2/)') 'Cluster #',i
 316           write (iout,'(3(1pe14.5),5x,1pe14.5)')
 317      &    ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
 318         endif
 319         W1i=0.0D0
 320         do k=2,3
 321           do l=2,3
 322             W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
 323           enddo
 324         enddo
 325         W1i=a(1,1)-W1i
 326         W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
 327 c        if (lprint) write(*,'(a,3(1pe15.5)/)')
 328 c     &          'detAp, W1, anormi',detApi,W1i,anormi
 329         do k=2,3
 330           zk=censc(k,i,it)
 331           do l=2,3
 332             zk=zk+zz1*Ap_inv(k,l)*a(l,1)
 333           enddo
 334           z(k,i)=zk
 335         enddo
 336         detAp(i)=dsqrt(detApi)
 337       enddo
 338
 339       if (lprint) then
 340         print *,'W1:',(w1(i),i=1,nlobit)
 341         print *,'detAp:',(detAp(i),i=1,nlobit)
 342         print *,'Z'
 343         do i=1,nlobit
 344           print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
 345         enddo
 346         write (iout,*) 'W1:',(w1(i),i=1,nlobit)
 347         write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
 348         write (iout,*) 'Z'
 349         do i=1,nlobit
 350           write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
 351         enddo
 352       endif
 353       if (lcheck) then
 354 C Writing the distribution just to check the procedure
 355       fac=0.0D0
 356       dV=deg2rad**2*10.0D0
 357       sum=0.0D0
 358       sum1=0.0D0
 359       do i=1,nlobit
 360         fac=fac+W1(i)/detAp(i)
 361       enddo
 362       fac=1.0D0/(2.0D0*fac*pi)
 363 cd    print *,it,'fac=',fac
 364       do ial=90,180,2
 365         y(1)=deg2rad*ial
 366         do iom=-180,180,5
 367           y(2)=deg2rad*iom
 368           wart=0.0D0
 369           do i=1,nlobit
 370             do j=2,3
 371               do k=2,3
 372                 a(j-1,k-1)=gaussc(j,k,i,it)
 373               enddo
 374             enddo
 375             y2=y(2)
 376
 377             do iii=-1,1
 378
 379               y(2)=y2+iii*dwapi
 380
 381               wykl=0.0D0
 382               do j=1,2
 383                 do k=1,2
 384                   wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
 385                 enddo
 386               enddo
 387               wart=wart+W1(i)*dexp(-0.5D0*wykl)
 388
 389             enddo
 390
 391             y(2)=y2
 392
 393           enddo
 394 c         print *,'y',y(1),y(2),' fac=',fac
 395           wart=fac*wart
 396           write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
 397           sum=sum+wart
 398           sum1=sum1+1.0D0
 399         enddo
 400       enddo
 401 c     print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
 402       return
 403       endif
 404
 405 C Calculate the CM of the system
 406 C
 407       do i=1,nlobit
 408         W1(i)=W1(i)/detAp(i)
 409       enddo
 410       sumW(0)=0.0D0
 411       do i=1,nlobit
 412         sumW(i)=sumW(i-1)+W1(i)
 413       enddo
 414       cm(1)=z(2,1)*W1(1)
 415       cm(2)=z(3,1)*W1(1)
 416       do j=2,nlobit
 417         cm(1)=cm(1)+z(2,j)*W1(j)
 418         cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
 419       enddo
 420       cm(1)=cm(1)/sumW(nlobit)
 421       cm(2)=cm(2)/sumW(nlobit)
 422       if (cm(1).gt.Big .or. cm(1).lt.-Big .or.
 423      & cm(2).gt.Big .or. cm(2).lt.-Big) then
 424 cd        write (iout,'(a)')
 425 cd     & 'Unexpected error in GenSide - CM coordinates too large.'
 426 cd        write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
 427 cd        write (*,'(a)')
 428 cd     & 'Unexpected error in GenSide - CM coordinates too large.'
 429 cd        write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
 430         fail=.true.
 431         return
 432       endif
 433 cd    print *,'CM:',cm(1),cm(2)
 434 C
 435 C Find the largest search distance from CM
 436 C
 437       radmax=0.0D0
 438       do i=1,nlobit
 439         do j=2,3
 440           do k=2,3
 441             a(j-1,k-1)=gaussc(j,k,i,it)
 442           enddo
 443         enddo
 444 #ifdef NAG
 445         call f02faf('N','U',2,a,3,eig,work,100,ifail)
 446 #else
 447         call djacob(2,3,10000,1.0d-10,a,vec,eig)
 448 #endif
 449 #ifdef MPI
 450         if (lprint) then
 451           print *,'*************** CG Processor',me
 452           print *,'CM:',cm(1),cm(2)
 453           write (iout,*) '*************** CG Processor',me
 454           write (iout,*) 'CM:',cm(1),cm(2)
 455           print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
 456           write (iout,'(A,8f10.5)')
 457      &        'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
 458         endif
 459 #endif
 460         if (eig(1).lt.eig_limit) then
 461           write(iout,'(a)')
 462      &     'From Mult_Norm: Eigenvalues of A are too small.'
 463           write(*,'(a)')
 464      &     'From Mult_Norm: Eigenvalues of A are too small.'
 465           fail=.true.
 466           return
 467         endif
 468         radius=0.0D0
 469 cd      print *,'i=',i
 470         do j=1,2
 471           radius=radius+pinorm(z(j+1,i)-cm(j))**2
 472         enddo
 473         radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
 474         if (radius.gt.radmax) radmax=radius
 475       enddo
 476       if (radmax.gt.pi) radmax=pi
 477 C
 478 C Determine the boundaries of the search rectangle.
 479 C
 480       if (lprint) then
 481         print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
 482         print '(a,4(1pe14.4))','radmax: ',radmax
 483       endif
 484       box(1,1)=dmax1(cm(1)-radmax,0.0D0)
 485       box(2,1)=dmin1(cm(1)+radmax,pi)
 486       box(1,2)=cm(2)-radmax
 487       box(2,2)=cm(2)+radmax
 488       if (lprint) then
 489 #ifdef MPI
 490         print *,'CG Processor',me,' Array BOX:'
 491 #else
 492         print *,'Array BOX:'
 493 #endif
 494         print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
 495         print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
 496 #ifdef MPI
 497         write (iout,*)'CG Processor',me,' Array BOX:'
 498 #else
 499         write (iout,*)'Array BOX:'
 500 #endif
 501         write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
 502         write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
 503       endif
 504       if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
 505 #ifdef MPI
 506         write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
 507         write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
 508 #else
 509 c        write (iout,'(a)') 'Bad sampling box.'
 510 #endif
 511         fail=.true.
 512         return
 513       endif
 514       which_lobe=ran_number(0.0D0,sumW(nlobit))
 515 c     print '(a,1pe14.4)','which_lobe=',which_lobe
 516       do i=1,nlobit
 517         if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
 518       enddo
 519     1 ilob=i
 520 c     print *,'ilob=',ilob,' nlob=',nlob(it)
 521       do i=2,3
 522         cm(i-1)=z(i,ilob)
 523         do j=2,3
 524           a(i-1,j-1)=gaussc(i,j,ilob,it)
 525         enddo
 526       enddo
 527 cd    print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
 528       call mult_norm1(3,2,a,cm,box,y,fail)
 529       if (fail) return
 530       al=y(1)
 531       om=pinorm(y(2))
 532 cd    print *,'al=',al,' om=',om
 533 cd    stop
 534       return
 535       end
 536 c---------------------------------------------------------------------------
 537       double precision function ran_number(x1,x2)
 538 C Calculate a random real number from the range (x1,x2).
 539       implicit real*8 (a-h,o-z)
 540       include 'DIMENSIONS'
 541       double precision x1,x2,fctor
 542       data fctor /2147483647.0D0/
 543 #ifdef MPI
 544       include "mpif.h"
 545       include 'COMMON.SETUP'
 546       ran_number=x1+(x2-x1)*prng_next(me)
 547 #else
 548       call vrnd(ix,1)
 549       ran_number=x1+(x2-x1)*ix/fctor
 550 #endif
 551       return
 552       end
 553 c--------------------------------------------------------------------------
 554       integer function iran_num(n1,n2)
 555 C Calculate a random integer number from the range (n1,n2).
 556       implicit real*8 (a-h,o-z)
 557       include 'DIMENSIONS'
 558       integer n1,n2,ix
 559       real fctor /2147483647.0/
 560 #ifdef MPI
 561       include "mpif.h"
 562       include 'COMMON.SETUP'
 563       ix=n1+(n2-n1+1)*prng_next(me)
 564       if (ix.lt.n1) ix=n1
 565       if (ix.gt.n2) ix=n2
 566       iran_num=ix
 567 #else
 568       call vrnd(ix,1)
 569       ix=n1+(n2-n1+1)*(ix/fctor)
 570       if (ix.gt.n2) ix=n2
 571       iran_num=ix
 572 #endif
 573       return
 574       end
 575 c--------------------------------------------------------------------------
 576       double precision function binorm(x1,x2,sigma1,sigma2,ak)
 577       implicit real*8 (a-h,o-z)
 578 c     print '(a)','Enter BINORM.'
 579       alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
 580       aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
 581       seg=sigma1/(sigma1+ak*sigma2)
 582       alen=ran_number(0.0D0,1.0D0)
 583       if (alen.lt.seg) then
 584         binorm=anorm_distr(x1,sigma1,alowb,aupb)
 585       else
 586         binorm=anorm_distr(x2,sigma2,alowb,aupb)
 587       endif
 588 c     print '(a)','Exiting BINORM.'
 589       return
 590       end
 591 c-----------------------------------------------------------------------
 592 c      double precision function anorm_distr(x,sigma,alowb,aupb)
 593 c      implicit real*8 (a-h,o-z)
 594 c     print '(a)','Enter ANORM_DISTR.'
 595 c   10 y=ran_number(alowb,aupb)
 596 c      expon=dexp(-0.5D0*((y-x)/sigma)**2)
 597 c      ran=ran_number(0.0D0,1.0D0)
 598 c      if (expon.lt.ran) goto 10
 599 c      anorm_distr=y
 600 c     print '(a)','Exiting ANORM_DISTR.'
 601 c      return
 602 c      end
 603 c-----------------------------------------------------------------------
 604         double precision function anorm_distr(x,sigma,alowb,aupb)
 605         implicit real*8 (a-h,o-z)
 606 c  to make a normally distributed deviate with zero mean and unit variance
 607 c
 608         integer iset
 609         real fac,gset,rsq,v1,v2,ran1
 610         save iset,gset
 611         data iset/0/
 612         if(iset.eq.0) then
 613 1               v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
 614                 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
 615                 rsq=v1**2+v2**2
 616                 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
 617                 fac=sqrt(-2.0d0*log(rsq)/rsq)
 618                 gset=v1*fac
 619                 gaussdev=v2*fac
 620                 iset=1
 621         else
 622                 gaussdev=gset
 623                 iset=0
 624         endif
 625         anorm_distr=x+gaussdev*sigma
 626       return
 627       end
 628 c------------------------------------------------------------------------
 629       subroutine mult_norm(lda,n,a,x,fail)
 630 C
 631 C Generate the vector X whose elements obey the multiple-normal distribution
 632 C from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
 633 C n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
 634 C eigenvalue of the matrix A is close to 0.
 635 C
 636       implicit double precision (a-h,o-z)
 637       double precision a(lda,n),x(n),eig(100),vec(3,3),work(100)
 638       double precision eig_limit /1.0D-8/
 639       logical fail
 640       fail=.false.
 641 c     print '(a)','Enter MULT_NORM.'
 642 C
 643 C Find the smallest eigenvalue of the matrix A.
 644 C
 645 c     do i=1,n
 646 c       print '(8f10.5)',(a(i,j),j=1,n)
 647 c     enddo
 648 #ifdef NAG
 649       call f02faf('V','U',2,a,lda,eig,work,100,ifail)
 650 #else
 651       call djacob(2,lda,10000,1.0d-10,a,vec,eig)
 652 #endif
 653 c     print '(8f10.5)',(eig(i),i=1,n)
 654 C     print '(a)'
 655 c     do i=1,n
 656 c       print '(8f10.5)',(a(i,j),j=1,n)
 657 c     enddo
 658       if (eig(1).lt.eig_limit) then
 659         print *,'From Mult_Norm: Eigenvalues of A are too small.'
 660         fail=.true.
 661         return
 662       endif
 663 C
 664 C Generate points following the normal distributions along the principal
 665 C axes of the moment matrix. Store in WORK.
 666 C
 667       do i=1,n
 668         sigma=1.0D0/dsqrt(eig(i))
 669         alim=-3.0D0*sigma
 670         work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
 671       enddo
 672 C
 673 C Transform the vector of normal variables back to the original basis.
 674 C
 675       do i=1,n
 676         xi=0.0D0
 677         do j=1,n
 678           xi=xi+a(i,j)*work(j)
 679         enddo
 680         x(i)=xi
 681       enddo
 682       return
 683       end
 684 c------------------------------------------------------------------------
 685       subroutine mult_norm1(lda,n,a,z,box,x,fail)
 686 C
 687 C Generate the vector X whose elements obey the multi-gaussian multi-dimensional
 688 C distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
 689 C leading dimension of the moment matrix A, n is the dimension of the
 690 C distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
 691 C smallest eigenvalue of the matrix A is close to 0.
 692 C
 693       implicit real*8 (a-h,o-z)
 694       include 'DIMENSIONS'
 695 #ifdef MPI
 696       include 'mpif.h'
 697 #endif
 698       double precision a(lda,n),z(n),x(n),box(n,n)
 699       double precision etmp
 700       include 'COMMON.IOUNITS'
 701 #ifdef MP
 702       include 'COMMON.SETUP'
 703 #endif
 704       logical fail
 705 C
 706 C Generate points following the normal distributions along the principal
 707 C axes of the moment matrix. Store in WORK.
 708 C
 709 cd    print *,'CG Processor',me,' entered MultNorm1.'
 710 cd    print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
 711 cd    do i=1,n
 712 cd      print *,i,box(1,i),box(2,i)
 713 cd    enddo
 714       istep = 0
 715    10 istep = istep + 1
 716       if (istep.gt.10000) then
 717 c        write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
 718 c     & ' in MultNorm1.'
 719 c        write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
 720 c     & ' in MultNorm1.'
 721 c        write (iout,*) 'box',box
 722 c        write (iout,*) 'a',a
 723 c        write (iout,*) 'z',z
 724         fail=.true.
 725         return
 726       endif
 727       do i=1,n
 728        x(i)=ran_number(box(1,i),box(2,i))
 729       enddo
 730       ww=0.0D0
 731       do i=1,n
 732         xi=pinorm(x(i)-z(i))
 733         ww=ww+0.5D0*a(i,i)*xi*xi
 734         do j=i+1,n
 735           ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
 736         enddo
 737       enddo
 738       dec=ran_number(0.0D0,1.0D0)
 739 c      print *,(x(i),i=1,n),ww,dexp(-ww),dec
 740 crc   if (dec.gt.dexp(-ww)) goto 10
 741       if(-ww.lt.100) then
 742        etmp=dexp(-ww)
 743       else
 744        return
 745       endif
 746       if (dec.gt.etmp) goto 10
 747 cd    print *,'CG Processor',me,' exitting MultNorm1.'
 748       return
 749       end
 750 c
 751 crc--------------------------------------
 752       subroutine overlap_sc(scfail)
 753 c     Internal and cartesian coordinates must be consistent as input,
 754 c     and will be up-to-date on return.
 755 c     At the end of this procedure, scfail is true if there are
 756 c     overlapping residues left, or false otherwise (success)
 757       implicit real*8 (a-h,o-z)
 758       include 'DIMENSIONS'
 759       include 'COMMON.CHAIN'
 760       include 'COMMON.INTERACT'
 761       include 'COMMON.FFIELD'
 762       include 'COMMON.VAR'
 763       include 'COMMON.SBRIDGE'
 764       include 'COMMON.IOUNITS'
 765       logical had_overlaps,fail,scfail
 766       integer ioverlap(maxres),ioverlap_last
 767
 768       had_overlaps=.false.
 769       call overlap_sc_list(ioverlap,ioverlap_last)
 770       if (ioverlap_last.gt.0) then
 771         write (iout,*) '#OVERLAPing residues ',ioverlap_last
 772         write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
 773         had_overlaps=.true.
 774       endif
 775
 776       maxsi=1000
 777       do k=1,1000
 778         if (ioverlap_last.eq.0) exit
 779
 780         do ires=1,ioverlap_last
 781           i=ioverlap(ires)
 782           iti=itype(i)
 783           if (iti.ne.10) then
 784             nsi=0
 785             fail=.true.
 786             do while (fail.and.nsi.le.maxsi)
 787               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 788               nsi=nsi+1
 789             enddo
 790             if(fail) goto 999
 791           endif
 792         enddo
 793
 794         call chainbuild
 795         call overlap_sc_list(ioverlap,ioverlap_last)
 796 c        write (iout,*) 'Overlaping residues ',ioverlap_last,
 797 c     &           (ioverlap(j),j=1,ioverlap_last)
 798       enddo
 799
 800       if (k.le.1000.and.ioverlap_last.eq.0) then
 801         scfail=.false.
 802         if (had_overlaps) then
 803           write (iout,*) '#OVERLAPing all corrected after ',k,
 804      &         ' random generation'
 805         endif
 806       else
 807         scfail=.true.
 808         write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
 809         write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
 810       endif
 811
 812       return
 813
 814  999  continue
 815       write (iout,'(a30,i5,a12,i4)')
 816      &               '#OVERLAP FAIL in gen_side after',maxsi,
 817      &               'iter for RES',i
 818       scfail=.true.
 819       return
 820       end
 821
 822       subroutine overlap_sc_list(ioverlap,ioverlap_last)
 823       implicit real*8 (a-h,o-z)
 824       include 'DIMENSIONS'
 825       include 'COMMON.GEO'
 826       include 'COMMON.LOCAL'
 827       include 'COMMON.IOUNITS'
 828       include 'COMMON.CHAIN'
 829       include 'COMMON.INTERACT'
 830       include 'COMMON.FFIELD'
 831       include 'COMMON.VAR'
 832       include 'COMMON.CALC'
 833       logical fail
 834       integer ioverlap(maxres),ioverlap_last
 835       data redfac /0.5D0/
 836
 837       ioverlap_last=0
 838 C Check for SC-SC overlaps and mark residues
 839 c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
 840       ind=0
 841       do i=iatsc_s,iatsc_e
 842         itypi=itype(i)
 843         itypi1=itype(i+1)
 844         xi=c(1,nres+i)
 845         yi=c(2,nres+i)
 846         zi=c(3,nres+i)
 847         dxi=dc_norm(1,nres+i)
 848         dyi=dc_norm(2,nres+i)
 849         dzi=dc_norm(3,nres+i)
 850         dsci_inv=dsc_inv(itypi)
 851 c
 852        do iint=1,nint_gr(i)
 853          do j=istart(i,iint),iend(i,iint)
 854             ind=ind+1
 855             itypj=itype(j)
 856             dscj_inv=dsc_inv(itypj)
 857             sig0ij=sigma(itypi,itypj)
 858             chi1=chi(itypi,itypj)
 859             chi2=chi(itypj,itypi)
 860             chi12=chi1*chi2
 861             chip1=chip(itypi)
 862             chip2=chip(itypj)
 863             chip12=chip1*chip2
 864             alf1=alp(itypi)
 865             alf2=alp(itypj)
 866             alf12=0.5D0*(alf1+alf2)
 867           if (j.gt.i+1) then
 868            rcomp=sigmaii(itypi,itypj)
 869           else
 870            rcomp=sigma(itypi,itypj)
 871           endif
 872 c         print '(2(a3,2i3),a3,2f10.5)',
 873 c     &        ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
 874 c     &        ,rcomp
 875             xj=c(1,nres+j)-xi
 876             yj=c(2,nres+j)-yi
 877             zj=c(3,nres+j)-zi
 878             dxj=dc_norm(1,nres+j)
 879             dyj=dc_norm(2,nres+j)
 880             dzj=dc_norm(3,nres+j)
 881             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 882             rij=dsqrt(rrij)
 883             call sc_angular
 884             sigsq=1.0D0/sigsq
 885             sig=sig0ij*dsqrt(sigsq)
 886             rij_shift=1.0D0/rij-sig+sig0ij
 887
 888 ct          if ( 1.0/rij .lt. redfac*rcomp .or.
 889 ct     &       rij_shift.le.0.0D0 ) then
 890             if ( rij_shift.le.0.0D0 ) then
 891 cd           write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
 892 cd     &     'overlap SC-SC: i=',i,' j=',j,
 893 cd     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
 894 cd     &     rcomp,1.0/rij,rij_shift
 895           ioverlap_last=ioverlap_last+1
 896           ioverlap(ioverlap_last)=i
 897           do k=1,ioverlap_last-1
 898            if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
 899           enddo
 900           ioverlap_last=ioverlap_last+1
 901           ioverlap(ioverlap_last)=j
 902           do k=1,ioverlap_last-1
 903            if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
 904           enddo
 905          endif
 906         enddo
 907        enddo
 908       enddo
 909       return
 910       end