source/unres/src_MD-M-SAXS/gen_rand_conf.F

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