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

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