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

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