source/unres/src_MD-M/gen_rand_conf.F

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