source/unres/src_CSA_DiL/gen_rand_conf.F

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