source/wham/src-new/ssMD.F

   1 c----------------------------------------------------------------------------
   2       subroutine check_energies
   3 c      implicit none
   4
   5 c     Includes
   6       include 'DIMENSIONS'
   7       include 'COMMON.CHAIN'
   8       include 'COMMON.VAR'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.SBRIDGE'
  11       include 'COMMON.LOCAL'
  12       include 'COMMON.GEO'
  13
  14 c     External functions
  15       double precision ran_number
  16       external ran_number
  17
  18 c     Local variables
  19       integer i,j,k,l,lmax,p,pmax
  20       double precision rmin,rmax
  21       double precision eij
  22
  23       double precision d
  24       double precision wi,rij,tj,pj
  25
  26
  27 c      return
  28
  29       i=5
  30       j=14
  31
  32       d=dsc(1)
  33       rmin=2.0D0
  34       rmax=12.0D0
  35
  36       lmax=10000
  37       pmax=1
  38
  39       do k=1,3
  40         c(k,i)=0.0D0
  41         c(k,j)=0.0D0
  42         c(k,nres+i)=0.0D0
  43         c(k,nres+j)=0.0D0
  44       enddo
  45
  46       do l=1,lmax
  47
  48 ct        wi=ran_number(0.0D0,pi)
  49 c        wi=ran_number(0.0D0,pi/6.0D0)
  50 c        wi=0.0D0
  51 ct        tj=ran_number(0.0D0,pi)
  52 ct        pj=ran_number(0.0D0,pi)
  53 c        pj=ran_number(0.0D0,pi/6.0D0)
  54 c        pj=0.0D0
  55
  56         do p=1,pmax
  57 ct           rij=ran_number(rmin,rmax)
  58
  59            c(1,j)=d*sin(pj)*cos(tj)
  60            c(2,j)=d*sin(pj)*sin(tj)
  61            c(3,j)=d*cos(pj)
  62
  63            c(3,nres+i)=-rij
  64
  65            c(1,i)=d*sin(wi)
  66            c(3,i)=-rij-d*cos(wi)
  67
  68            do k=1,3
  69               dc(k,nres+i)=c(k,nres+i)-c(k,i)
  70               dc_norm(k,nres+i)=dc(k,nres+i)/d
  71               dc(k,nres+j)=c(k,nres+j)-c(k,j)
  72               dc_norm(k,nres+j)=dc(k,nres+j)/d
  73            enddo
  74
  75            call dyn_ssbond_ene(i,j,eij)
  76         enddo
  77       enddo
  78
  79       call exit(1)
  80
  81       return
  82       end
  83
  84 C-----------------------------------------------------------------------------
  85
  86       subroutine dyn_ssbond_ene(resi,resj,eij)
  87 c      implicit none
  88
  89 c     Includes
  90       include 'DIMENSIONS'
  91       include 'COMMON.SBRIDGE'
  92       include 'COMMON.CHAIN'
  93       include 'COMMON.DERIV'
  94       include 'COMMON.LOCAL'
  95       include 'COMMON.INTERACT'
  96       include 'COMMON.VAR'
  97       include 'COMMON.IOUNITS'
  98       include 'COMMON.CALC'
  99 #ifndef CLUST
 100 #ifndef WHAM
 101       include 'COMMON.MD'
 102 #endif
 103 #endif
 104
 105 c     External functions
 106       double precision h_base
 107       external h_base
 108
 109 c     Input arguments
 110       integer resi,resj
 111
 112 c     Output arguments
 113       double precision eij
 114
 115 c     Local variables
 116       logical havebond
 117 c      integer itypi,itypj,k,l
 118       double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
 119       double precision sig0ij,ljd,sig,fac,e1,e2
 120       double precision dcosom1(3),dcosom2(3),ed
 121       double precision pom1,pom2
 122       double precision ljA,ljB,ljXs
 123       double precision d_ljB(1:3)
 124       double precision ssA,ssB,ssC,ssXs
 125       double precision ssxm,ljxm,ssm,ljm
 126       double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
 127       double precision f1,f2,h1,h2,hd1,hd2
 128       double precision omega,delta_inv,deltasq_inv,fac1,fac2
 129 c-------FIRST METHOD
 130       double precision xm,d_xm(1:3)
 131 c-------END FIRST METHOD
 132 c-------SECOND METHOD
 133 c$$$      double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
 134 c-------END SECOND METHOD
 135
 136 c-------TESTING CODE
 137       logical checkstop,transgrad
 138       common /sschecks/ checkstop,transgrad
 139
 140       integer icheck,nicheck,jcheck,njcheck
 141       double precision echeck(-1:1),deps,ssx0,ljx0
 142 c-------END TESTING CODE
 143
 144
 145       i=resi
 146       j=resj
 147
 148       itypi=itype(i)
 149       dxi=dc_norm(1,nres+i)
 150       dyi=dc_norm(2,nres+i)
 151       dzi=dc_norm(3,nres+i)
 152       dsci_inv=vbld_inv(i+nres)
 153
 154       itypj=itype(j)
 155       xj=c(1,nres+j)-c(1,nres+i)
 156       yj=c(2,nres+j)-c(2,nres+i)
 157       zj=c(3,nres+j)-c(3,nres+i)
 158       dxj=dc_norm(1,nres+j)
 159       dyj=dc_norm(2,nres+j)
 160       dzj=dc_norm(3,nres+j)
 161       dscj_inv=vbld_inv(j+nres)
 162
 163       chi1=chi(itypi,itypj)
 164       chi2=chi(itypj,itypi)
 165       chi12=chi1*chi2
 166       chip1=chip(itypi)
 167       chip2=chip(itypj)
 168       chip12=chip1*chip2
 169       alf1=alp(itypi)
 170       alf2=alp(itypj)
 171       alf12=0.5D0*(alf1+alf2)
 172
 173       rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 174       rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
 175 c     The following are set in sc_angular
 176 c      erij(1)=xj*rij
 177 c      erij(2)=yj*rij
 178 c      erij(3)=zj*rij
 179 c      om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
 180 c      om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
 181 c      om12=dxi*dxj+dyi*dyj+dzi*dzj
 182       call sc_angular
 183       rij=1.0D0/rij  ! Reset this so it makes sense
 184
 185       sig0ij=sigma(itypi,itypj)
 186       sig=sig0ij*dsqrt(1.0D0/sigsq)
 187
 188       ljXs=sig-sig0ij
 189       ljA=eps1*eps2rt**2*eps3rt**2
 190       ljB=ljA*bb(itypi,itypj)
 191       ljA=ljA*aa(itypi,itypj)
 192       ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 193
 194       ssXs=d0cm
 195       deltat1=1.0d0-om1
 196       deltat2=1.0d0+om2
 197       deltat12=om2-om1+2.0d0
 198       cosphi=om12-om1*om2
 199       ssA=akcm
 200       ssB=akct*deltat12
 201       ssC=ss_depth
 202      &     +akth*(deltat1*deltat1+deltat2*deltat2)
 203      &     +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
 204       ssxm=ssXs-0.5D0*ssB/ssA
 205
 206 c-------TESTING CODE
 207 c$$$c     Some extra output
 208 c$$$      ssm=ssC-0.25D0*ssB*ssB/ssA
 209 c$$$      ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
 210 c$$$      ssx0=ssB*ssB-4.0d0*ssA*ssC
 211 c$$$      if (ssx0.gt.0.0d0) then
 212 c$$$        ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
 213 c$$$      else
 214 c$$$        ssx0=ssxm
 215 c$$$      endif
 216 c$$$      ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 217 c$$$      write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
 218 c$$$     &     ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
 219 c$$$      return
 220 c-------END TESTING CODE
 221
 222 c-------TESTING CODE
 223 c     Stop and plot energy and derivative as a function of distance
 224       if (checkstop) then
 225         ssm=ssC-0.25D0*ssB*ssB/ssA
 226         ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
 227         if (ssm.lt.ljm .and.
 228      &       dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
 229           nicheck=1000
 230           njcheck=1
 231           deps=0.5d-7
 232         else
 233           checkstop=.false.
 234         endif
 235       endif
 236       if (.not.checkstop) then
 237         nicheck=0
 238         njcheck=-1
 239       endif
 240
 241       do icheck=0,nicheck
 242       do jcheck=-1,njcheck
 243       if (checkstop) rij=(ssxm-1.0d0)+
 244      &       ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
 245 c-------END TESTING CODE
 246
 247       if (rij.gt.ljxm) then
 248         havebond=.false.
 249         ljd=rij-ljXs
 250         fac=(1.0D0/ljd)**expon
 251         e1=fac*fac*aa(itypi,itypj)
 252         e2=fac*bb(itypi,itypj)
 253         eij=eps1*eps2rt*eps3rt*(e1+e2)
 254         eps2der=eij*eps3rt
 255         eps3der=eij*eps2rt
 256         eij=eij*eps2rt*eps3rt
 257
 258         sigder=-sig/sigsq
 259         e1=e1*eps1*eps2rt**2*eps3rt**2
 260         ed=-expon*(e1+eij)/ljd
 261         sigder=ed*sigder
 262         eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
 263         eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
 264         eom12=eij*eps1_om12+eps2der*eps2rt_om12
 265      &       -2.0D0*alf12*eps3der+sigder*sigsq_om12
 266       else if (rij.lt.ssxm) then
 267         havebond=.true.
 268         ssd=rij-ssXs
 269         eij=ssA*ssd*ssd+ssB*ssd+ssC
 270
 271         ed=2*akcm*ssd+akct*deltat12
 272         pom1=akct*ssd
 273         pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
 274         eom1=-2*akth*deltat1-pom1-om2*pom2
 275         eom2= 2*akth*deltat2+pom1-om1*pom2
 276         eom12=pom2
 277       else
 278         omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
 279
 280         d_ssxm(1)=0.5D0*akct/ssA
 281         d_ssxm(2)=-d_ssxm(1)
 282         d_ssxm(3)=0.0D0
 283
 284         d_ljxm(1)=sig0ij/sqrt(sigsq**3)
 285         d_ljxm(2)=d_ljxm(1)*sigsq_om2
 286         d_ljxm(3)=d_ljxm(1)*sigsq_om12
 287         d_ljxm(1)=d_ljxm(1)*sigsq_om1
 288
 289 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
 290         xm=0.5d0*(ssxm+ljxm)
 291         do k=1,3
 292           d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
 293         enddo
 294         if (rij.lt.xm) then
 295           havebond=.true.
 296           ssm=ssC-0.25D0*ssB*ssB/ssA
 297           d_ssm(1)=0.5D0*akct*ssB/ssA
 298           d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
 299           d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
 300           d_ssm(3)=omega
 301           f1=(rij-xm)/(ssxm-xm)
 302           f2=(rij-ssxm)/(xm-ssxm)
 303           h1=h_base(f1,hd1)
 304           h2=h_base(f2,hd2)
 305           eij=ssm*h1+Ht*h2
 306           delta_inv=1.0d0/(xm-ssxm)
 307           deltasq_inv=delta_inv*delta_inv
 308           fac=ssm*hd1-Ht*hd2
 309           fac1=deltasq_inv*fac*(xm-rij)
 310           fac2=deltasq_inv*fac*(rij-ssxm)
 311           ed=delta_inv*(Ht*hd2-ssm*hd1)
 312           eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
 313           eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
 314           eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
 315         else
 316           havebond=.false.
 317           ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
 318           d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
 319           d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
 320           d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
 321      +         alf12/eps3rt)
 322           d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
 323           f1=(rij-ljxm)/(xm-ljxm)
 324           f2=(rij-xm)/(ljxm-xm)
 325           h1=h_base(f1,hd1)
 326           h2=h_base(f2,hd2)
 327           eij=Ht*h1+ljm*h2
 328           delta_inv=1.0d0/(ljxm-xm)
 329           deltasq_inv=delta_inv*delta_inv
 330           fac=Ht*hd1-ljm*hd2
 331           fac1=deltasq_inv*fac*(ljxm-rij)
 332           fac2=deltasq_inv*fac*(rij-xm)
 333           ed=delta_inv*(ljm*hd2-Ht*hd1)
 334           eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
 335           eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
 336           eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
 337         endif
 338 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
 339
 340 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 341 c$$$        ssd=rij-ssXs
 342 c$$$        ljd=rij-ljXs
 343 c$$$        fac1=rij-ljxm
 344 c$$$        fac2=rij-ssxm
 345 c$$$
 346 c$$$        d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
 347 c$$$        d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
 348 c$$$        d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
 349 c$$$
 350 c$$$        ssm=ssC-0.25D0*ssB*ssB/ssA
 351 c$$$        d_ssm(1)=0.5D0*akct*ssB/ssA
 352 c$$$        d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
 353 c$$$        d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
 354 c$$$        d_ssm(3)=omega
 355 c$$$
 356 c$$$        ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
 357 c$$$        do k=1,3
 358 c$$$          d_ljm(k)=ljm*d_ljB(k)
 359 c$$$        enddo
 360 c$$$        ljm=ljm*ljB
 361 c$$$
 362 c$$$        ss=ssA*ssd*ssd+ssB*ssd+ssC
 363 c$$$        d_ss(0)=2.0d0*ssA*ssd+ssB
 364 c$$$        d_ss(2)=akct*ssd
 365 c$$$        d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
 366 c$$$        d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
 367 c$$$        d_ss(3)=omega
 368 c$$$
 369 c$$$        ljf=bb(itypi,itypj)/aa(itypi,itypj)
 370 c$$$        ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
 371 c$$$        d_ljf(0)=ljf*2.0d0*ljB*fac1
 372 c$$$        do k=1,3
 373 c$$$          d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
 374 c$$$     &         2.0d0*ljB*fac1*d_ljxm(k))
 375 c$$$        enddo
 376 c$$$        ljf=ljm+ljf*ljB*fac1*fac1
 377 c$$$
 378 c$$$        f1=(rij-ljxm)/(ssxm-ljxm)
 379 c$$$        f2=(rij-ssxm)/(ljxm-ssxm)
 380 c$$$        h1=h_base(f1,hd1)
 381 c$$$        h2=h_base(f2,hd2)
 382 c$$$        eij=ss*h1+ljf*h2
 383 c$$$        delta_inv=1.0d0/(ljxm-ssxm)
 384 c$$$        deltasq_inv=delta_inv*delta_inv
 385 c$$$        fac=ljf*hd2-ss*hd1
 386 c$$$        ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
 387 c$$$        eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
 388 c$$$     &       (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
 389 c$$$        eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
 390 c$$$     &       (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
 391 c$$$        eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
 392 c$$$     &       (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
 393 c$$$
 394 c$$$        havebond=.false.
 395 c$$$        if (ed.gt.0.0d0) havebond=.true.
 396 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 397
 398       endif
 399
 400       if (havebond) then
 401 #ifndef CLUST
 402 #ifndef WHAM
 403 c        if (dyn_ssbond_ij(i,j).eq.1.0d300) then
 404 c          write(iout,'(a15,f12.2,f8.1,2i5)')
 405 c     &         "SSBOND_E_FORM",totT,t_bath,i,j
 406 c        endif
 407 #endif
 408 #endif
 409       write(iout,*), 'DYN_SS_BOND',i,j,eij
 410         dyn_ssbond_ij(i,j)=eij
 411       else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
 412         dyn_ssbond_ij(i,j)=1.0d300
 413 #ifndef CLUST
 414 #ifndef WHAM
 415 c        write(iout,'(a15,f12.2,f8.1,2i5)')
 416 c     &       "SSBOND_E_BREAK",totT,t_bath,i,j
 417 #endif
 418 #endif
 419       endif
 420
 421 c-------TESTING CODE
 422       if (checkstop) then
 423         if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
 424      &       "CHECKSTOP",rij,eij,ed
 425         echeck(jcheck)=eij
 426       endif
 427       enddo
 428       if (checkstop) then
 429         write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
 430       endif
 431       enddo
 432       if (checkstop) then
 433         transgrad=.true.
 434         checkstop=.false.
 435       endif
 436 c-------END TESTING CODE
 437
 438       do k=1,3
 439         dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
 440         dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
 441       enddo
 442       do k=1,3
 443         gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
 444       enddo
 445       do k=1,3
 446         gvdwx(k,i)=gvdwx(k,i)-gg(k)
 447      &       +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
 448      &       +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 449         gvdwx(k,j)=gvdwx(k,j)+gg(k)
 450      &       +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
 451      &       +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
 452       enddo
 453 cgrad      do k=i,j-1
 454 cgrad        do l=1,3
 455 cgrad          gvdwc(l,k)=gvdwc(l,k)+gg(l)
 456 cgrad        enddo
 457 cgrad      enddo
 458
 459       do l=1,3
 460         gvdwc(l,i)=gvdwc(l,i)-gg(l)
 461         gvdwc(l,j)=gvdwc(l,j)+gg(l)
 462       enddo
 463
 464       return
 465       end
 466
 467 C-----------------------------------------------------------------------------
 468
 469       double precision function h_base(x,deriv)
 470 c     A smooth function going 0->1 in range [0,1]
 471 c     It should NOT be called outside range [0,1], it will not work there.
 472       implicit none
 473
 474 c     Input arguments
 475       double precision x
 476
 477 c     Output arguments
 478       double precision deriv
 479
 480 c     Local variables
 481       double precision xsq
 482
 483
 484 c     Two parabolas put together.  First derivative zero at extrema
 485 c$$$      if (x.lt.0.5D0) then
 486 c$$$        h_base=2.0D0*x*x
 487 c$$$        deriv=4.0D0*x
 488 c$$$      else
 489 c$$$        deriv=1.0D0-x
 490 c$$$        h_base=1.0D0-2.0D0*deriv*deriv
 491 c$$$        deriv=4.0D0*deriv
 492 c$$$      endif
 493
 494 c     Third degree polynomial.  First derivative zero at extrema
 495       h_base=x*x*(3.0d0-2.0d0*x)
 496       deriv=6.0d0*x*(1.0d0-x)
 497
 498 c     Fifth degree polynomial.  First and second derivatives zero at extrema
 499 c$$$      xsq=x*x
 500 c$$$      h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
 501 c$$$      deriv=x-1.0d0
 502 c$$$      deriv=deriv*deriv
 503 c$$$      deriv=30.0d0*xsq*deriv
 504
 505       return
 506       end
 507
 508 c----------------------------------------------------------------------------
 509
 510       subroutine dyn_set_nss
 511 c     Adjust nss and other relevant variables based on dyn_ssbond_ij
 512 c      implicit none
 513
 514 c     Includes
 515       include 'DIMENSIONS'
 516 #ifdef MPI
 517       include "mpif.h"
 518 #endif
 519       include 'COMMON.SBRIDGE'
 520       include 'COMMON.CHAIN'
 521       include 'COMMON.IOUNITS'
 522 c#ifdef CLUST
 523 c#ifdef WHAM
 524       integer max_fg_procs,maxprocs
 525       parameter (max_fg_procs=128)
 526       parameter (maxprocs=128)
 527 c#endif
 528 c#endif
 529 cc      include 'COMMON.SETUP'
 530 #ifndef CLUST
 531 #ifndef WHAM
 532       include 'COMMON.MD'
 533       include 'COMMON.SETUP'
 534 #endif
 535 #endif
 536
 537 c     Local variables
 538       double precision emin
 539       integer i,j,imin
 540       integer diff,allflag(maxdim),allnss,
 541      &     allihpb(maxdim),alljhpb(maxdim),
 542      &     newnss,newihpb(maxdim),newjhpb(maxdim)
 543       logical found
 544       integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
 545       integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
 546
 547       allnss=0
 548       do i=1,nres-1
 549         do j=i+1,nres
 550           if (dyn_ssbond_ij(i,j).lt.1.0d300) then
 551             allnss=allnss+1
 552             allflag(allnss)=0
 553             allihpb(allnss)=i
 554             alljhpb(allnss)=j
 555           endif
 556         enddo
 557       enddo
 558
 559 cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
 560
 561  1    emin=1.0d300
 562       do i=1,allnss
 563         if (allflag(i).eq.0 .and.
 564      &       dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
 565           emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
 566           imin=i
 567         endif
 568       enddo
 569       if (emin.lt.1.0d300) then
 570         allflag(imin)=1
 571         do i=1,allnss
 572           if (allflag(i).eq.0 .and.
 573      &         (allihpb(i).eq.allihpb(imin) .or.
 574      &         alljhpb(i).eq.allihpb(imin) .or.
 575      &         allihpb(i).eq.alljhpb(imin) .or.
 576      &         alljhpb(i).eq.alljhpb(imin))) then
 577             allflag(i)=-1
 578           endif
 579         enddo
 580         goto 1
 581       endif
 582
 583 cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
 584
 585       newnss=0
 586       do i=1,allnss
 587         if (allflag(i).eq.1) then
 588           newnss=newnss+1
 589           newihpb(newnss)=allihpb(i)
 590           newjhpb(newnss)=alljhpb(i)
 591         endif
 592       enddo
 593
 594 #ifdef MPI
 595       if (nfgtasks.gt.1)then
 596
 597         call MPI_Reduce(newnss,g_newnss,1,
 598      &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
 599         call MPI_Gather(newnss,1,MPI_INTEGER,
 600      &                  i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
 601         displ(0)=0
 602         do i=1,nfgtasks-1,1
 603           displ(i)=i_newnss(i-1)+displ(i-1)
 604         enddo
 605         call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
 606      &                   g_newihpb,i_newnss,displ,MPI_INTEGER,
 607      &                   king,FG_COMM,IERR)
 608         call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
 609      &                   g_newjhpb,i_newnss,displ,MPI_INTEGER,
 610      &                   king,FG_COMM,IERR)
 611         if(fg_rank.eq.0) then
 612 c         print *,'g_newnss',g_newnss
 613 c         print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
 614 c         print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
 615          newnss=g_newnss
 616          do i=1,newnss
 617           newihpb(i)=g_newihpb(i)
 618           newjhpb(i)=g_newjhpb(i)
 619          enddo
 620         endif
 621       endif
 622 #endif
 623
 624       diff=newnss-nss
 625
 626 cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
 627
 628       do i=1,nss
 629         found=.false.
 630         do j=1,newnss
 631           if (idssb(i).eq.newihpb(j) .and.
 632      &         jdssb(i).eq.newjhpb(j)) found=.true.
 633         enddo
 634 #ifndef CLUST
 635 #ifndef WHAM
 636         if (.not.found.and.fg_rank.eq.0)
 637      &      write(iout,'(a15,f12.2,f8.1,2i5)')
 638      &       "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
 639 #endif
 640 #endif
 641       enddo
 642
 643       do i=1,newnss
 644         found=.false.
 645         do j=1,nss
 646           if (newihpb(i).eq.idssb(j) .and.
 647      &         newjhpb(i).eq.jdssb(j)) found=.true.
 648         enddo
 649 #ifndef CLUST
 650 #ifndef WHAM
 651         if (.not.found.and.fg_rank.eq.0)
 652      &      write(iout,'(a15,f12.2,f8.1,2i5)')
 653      &       "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
 654 #endif
 655 #endif
 656       enddo
 657
 658       nss=newnss
 659       do i=1,nss
 660         idssb(i)=newihpb(i)
 661         jdssb(i)=newjhpb(i)
 662       enddo
 663
 664       return
 665       end
 666
 667 c----------------------------------------------------------------------------
 668 #ifdef NIEWIADOM
 669 c#ifdef WHAM
 670       subroutine read_ssHist
 671       implicit none
 672
 673 c     Includes
 674       include 'DIMENSIONS'
 675       include "DIMENSIONS.FREE"
 676       include 'COMMON.FREE'
 677
 678 c     Local variables
 679       integer i,j
 680       character*80 controlcard
 681
 682       do i=1,dyn_nssHist
 683         call card_concat(controlcard,.true.)
 684         read(controlcard,*)
 685      &       dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
 686       enddo
 687
 688       return
 689       end
 690 #endif
 691
 692 c----------------------------------------------------------------------------
 693
 694
 695 C-----------------------------------------------------------------------------
 696 C-----------------------------------------------------------------------------
 697 C-----------------------------------------------------------------------------
 698 C-----------------------------------------------------------------------------
 699 C-----------------------------------------------------------------------------
 700 C-----------------------------------------------------------------------------
 701 C-----------------------------------------------------------------------------
 702
 703 c$$$c-----------------------------------------------------------------------------
 704 c$$$
 705 c$$$      subroutine ss_relax(i_in,j_in)
 706 c$$$      implicit none
 707 c$$$
 708 c$$$c     Includes
 709 c$$$      include 'DIMENSIONS'
 710 c$$$      include 'COMMON.VAR'
 711 c$$$      include 'COMMON.CHAIN'
 712 c$$$      include 'COMMON.IOUNITS'
 713 c$$$      include 'COMMON.INTERACT'
 714 c$$$
 715 c$$$c     Input arguments
 716 c$$$      integer i_in,j_in
 717 c$$$
 718 c$$$c     Local variables
 719 c$$$      integer i,iretcode,nfun_sc
 720 c$$$      logical scfail
 721 c$$$      double precision var(maxvar),e_sc,etot
 722 c$$$
 723 c$$$
 724 c$$$      mask_r=.true.
 725 c$$$      do i=nnt,nct
 726 c$$$        mask_side(i)=0
 727 c$$$      enddo
 728 c$$$      mask_side(i_in)=1
 729 c$$$      mask_side(j_in)=1
 730 c$$$
 731 c$$$c     Minimize the two selected side-chains
 732 c$$$      call overlap_sc(scfail)  ! Better not fail!
 733 c$$$      call minimize_sc(e_sc,var,iretcode,nfun_sc)
 734 c$$$
 735 c$$$      mask_r=.false.
 736 c$$$
 737 c$$$      return
 738 c$$$      end
 739 c$$$
 740 c$$$c-------------------------------------------------------------
 741 c$$$
 742 c$$$      subroutine minimize_sc(etot_sc,iretcode,nfun)
 743 c$$$c     Minimize side-chains only, starting from geom but without modifying
 744 c$$$c     bond lengths.
 745 c$$$c     If mask_r is already set, only the selected side-chains are minimized,
 746 c$$$c     otherwise all side-chains are minimized keeping the backbone frozen.
 747 c$$$      implicit none
 748 c$$$
 749 c$$$c     Includes
 750 c$$$      include 'DIMENSIONS'
 751 c$$$      include 'COMMON.IOUNITS'
 752 c$$$      include 'COMMON.VAR'
 753 c$$$      include 'COMMON.CHAIN'
 754 c$$$      include 'COMMON.GEO'
 755 c$$$      include 'COMMON.MINIM'
 756 c$$$      integer icall
 757 c$$$      common /srutu/ icall
 758 c$$$
 759 c$$$c     Output arguments
 760 c$$$      double precision etot_sc
 761 c$$$      integer iretcode,nfun
 762 c$$$
 763 c$$$c     External functions/subroutines
 764 c$$$      external func_sc,grad_sc,fdum
 765 c$$$
 766 c$$$c     Local variables
 767 c$$$      integer liv,lv
 768 c$$$      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 769 c$$$      integer iv(liv)
 770 c$$$      double precision rdum(1)
 771 c$$$      double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
 772 c$$$      integer idum(1)
 773 c$$$      integer i,nvar_restr
 774 c$$$
 775 c$$$
 776 c$$$cmc      start_minim=.true.
 777 c$$$      call deflt(2,iv,liv,lv,v)
 778 c$$$* 12 means fresh start, dont call deflt
 779 c$$$      iv(1)=12
 780 c$$$* max num of fun calls
 781 c$$$      if (maxfun.eq.0) maxfun=500
 782 c$$$      iv(17)=maxfun
 783 c$$$* max num of iterations
 784 c$$$      if (maxmin.eq.0) maxmin=1000
 785 c$$$      iv(18)=maxmin
 786 c$$$* controls output
 787 c$$$      iv(19)=1
 788 c$$$* selects output unit
 789 c$$$      iv(21)=0
 790 c$$$c      iv(21)=iout               ! DEBUG
 791 c$$$c      iv(21)=8                  ! DEBUG
 792 c$$$* 1 means to print out result
 793 c$$$      iv(22)=0
 794 c$$$c      iv(22)=1                  ! DEBUG
 795 c$$$* 1 means to print out summary stats
 796 c$$$      iv(23)=0
 797 c$$$c      iv(23)=1                  ! DEBUG
 798 c$$$* 1 means to print initial x and d
 799 c$$$      iv(24)=0
 800 c$$$c      iv(24)=1                  ! DEBUG
 801 c$$$* min val for v(radfac) default is 0.1
 802 c$$$      v(24)=0.1D0
 803 c$$$* max val for v(radfac) default is 4.0
 804 c$$$      v(25)=2.0D0
 805 c$$$c     v(25)=4.0D0
 806 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 807 c$$$* the sumsl default is 0.1
 808 c$$$      v(26)=0.1D0
 809 c$$$* false conv if (act fnctn decrease) .lt. v(34)
 810 c$$$* the sumsl default is 100*machep
 811 c$$$      v(34)=v(34)/100.0D0
 812 c$$$* absolute convergence
 813 c$$$      if (tolf.eq.0.0D0) tolf=1.0D-4
 814 c$$$      v(31)=tolf
 815 c$$$* relative convergence
 816 c$$$      if (rtolf.eq.0.0D0) rtolf=1.0D-1
 817 c$$$      v(32)=rtolf
 818 c$$$* controls initial step size
 819 c$$$       v(35)=1.0D-1
 820 c$$$* large vals of d correspond to small components of step
 821 c$$$      do i=1,nphi
 822 c$$$        d(i)=1.0D-1
 823 c$$$      enddo
 824 c$$$      do i=nphi+1,nvar
 825 c$$$        d(i)=1.0D-1
 826 c$$$      enddo
 827 c$$$
 828 c$$$      call geom_to_var(nvar,x)
 829 c$$$      IF (mask_r) THEN
 830 c$$$        do i=1,nres             ! Just in case...
 831 c$$$          mask_phi(i)=0
 832 c$$$          mask_theta(i)=0
 833 c$$$        enddo
 834 c$$$        call x2xx(x,xx,nvar_restr)
 835 c$$$        call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
 836 c$$$     &       iv,liv,lv,v,idum,rdum,fdum)
 837 c$$$        call xx2x(x,xx)
 838 c$$$      ELSE
 839 c$$$c     When minimizing ALL side-chains, etotal_sc is a little
 840 c$$$c     faster if we don't set mask_r
 841 c$$$        do i=1,nres
 842 c$$$          mask_phi(i)=0
 843 c$$$          mask_theta(i)=0
 844 c$$$          mask_side(i)=1
 845 c$$$        enddo
 846 c$$$        call x2xx(x,xx,nvar_restr)
 847 c$$$        call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
 848 c$$$     &       iv,liv,lv,v,idum,rdum,fdum)
 849 c$$$        call xx2x(x,xx)
 850 c$$$      ENDIF
 851 c$$$      call var_to_geom(nvar,x)
 852 c$$$      call chainbuild_sc
 853 c$$$      etot_sc=v(10)
 854 c$$$      iretcode=iv(1)
 855 c$$$      nfun=iv(6)
 856 c$$$      return
 857 c$$$      end
 858 c$$$
 859 c$$$C--------------------------------------------------------------------------
 860 c$$$
 861 c$$$      subroutine chainbuild_sc
 862 c$$$      implicit none
 863 c$$$      include 'DIMENSIONS'
 864 c$$$      include 'COMMON.VAR'
 865 c$$$      include 'COMMON.INTERACT'
 866 c$$$
 867 c$$$c     Local variables
 868 c$$$      integer i
 869 c$$$
 870 c$$$
 871 c$$$      do i=nnt,nct
 872 c$$$        if (.not.mask_r .or. mask_side(i).eq.1) then
 873 c$$$          call locate_side_chain(i)
 874 c$$$        endif
 875 c$$$      enddo
 876 c$$$
 877 c$$$      return
 878 c$$$      end
 879 c$$$
 880 c$$$C--------------------------------------------------------------------------
 881 c$$$
 882 c$$$      subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
 883 c$$$      implicit none
 884 c$$$
 885 c$$$c     Includes
 886 c$$$      include 'DIMENSIONS'
 887 c$$$      include 'COMMON.DERIV'
 888 c$$$      include 'COMMON.VAR'
 889 c$$$      include 'COMMON.MINIM'
 890 c$$$      include 'COMMON.IOUNITS'
 891 c$$$
 892 c$$$c     Input arguments
 893 c$$$      integer n
 894 c$$$      double precision x(maxvar)
 895 c$$$      double precision ufparm
 896 c$$$      external ufparm
 897 c$$$
 898 c$$$c     Input/Output arguments
 899 c$$$      integer nf
 900 c$$$      integer uiparm(1)
 901 c$$$      double precision urparm(1)
 902 c$$$
 903 c$$$c     Output arguments
 904 c$$$      double precision f
 905 c$$$
 906 c$$$c     Local variables
 907 c$$$      double precision energia(0:n_ene)
 908 c$$$#ifdef OSF
 909 c$$$c     Variables used to intercept NaNs
 910 c$$$      double precision x_sum
 911 c$$$      integer i_NAN
 912 c$$$#endif
 913 c$$$
 914 c$$$
 915 c$$$      nfl=nf
 916 c$$$      icg=mod(nf,2)+1
 917 c$$$
 918 c$$$#ifdef OSF
 919 c$$$c     Intercept NaNs in the coordinates, before calling etotal_sc
 920 c$$$      x_sum=0.D0
 921 c$$$      do i_NAN=1,n
 922 c$$$        x_sum=x_sum+x(i_NAN)
 923 c$$$      enddo
 924 c$$$c     Calculate the energy only if the coordinates are ok
 925 c$$$      if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
 926 c$$$        write(iout,*)"   *** func_restr_sc : Found NaN in coordinates"
 927 c$$$        f=1.0D+77
 928 c$$$        nf=0
 929 c$$$      else
 930 c$$$#endif
 931 c$$$
 932 c$$$      call var_to_geom_restr(n,x)
 933 c$$$      call zerograd
 934 c$$$      call chainbuild_sc
 935 c$$$      call etotal_sc(energia(0))
 936 c$$$      f=energia(0)
 937 c$$$      if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
 938 c$$$
 939 c$$$#ifdef OSF
 940 c$$$      endif
 941 c$$$#endif
 942 c$$$
 943 c$$$      return
 944 c$$$      end
 945 c$$$
 946 c$$$c-------------------------------------------------------
 947 c$$$
 948 c$$$      subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
 949 c$$$      implicit none
 950 c$$$
 951 c$$$c     Includes
 952 c$$$      include 'DIMENSIONS'
 953 c$$$      include 'COMMON.CHAIN'
 954 c$$$      include 'COMMON.DERIV'
 955 c$$$      include 'COMMON.VAR'
 956 c$$$      include 'COMMON.INTERACT'
 957 c$$$      include 'COMMON.MINIM'
 958 c$$$
 959 c$$$c     Input arguments
 960 c$$$      integer n
 961 c$$$      double precision x(maxvar)
 962 c$$$      double precision ufparm
 963 c$$$      external ufparm
 964 c$$$
 965 c$$$c     Input/Output arguments
 966 c$$$      integer nf
 967 c$$$      integer uiparm(1)
 968 c$$$      double precision urparm(1)
 969 c$$$
 970 c$$$c     Output arguments
 971 c$$$      double precision g(maxvar)
 972 c$$$
 973 c$$$c     Local variables
 974 c$$$      double precision f,gphii,gthetai,galphai,gomegai
 975 c$$$      integer ig,ind,i,j,k,igall,ij
 976 c$$$
 977 c$$$
 978 c$$$      icg=mod(nf,2)+1
 979 c$$$      if (nf-nfl+1) 20,30,40
 980 c$$$   20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
 981 c$$$c     write (iout,*) 'grad 20'
 982 c$$$      if (nf.eq.0) return
 983 c$$$      goto 40
 984 c$$$   30 call var_to_geom_restr(n,x)
 985 c$$$      call chainbuild_sc
 986 c$$$C
 987 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 988 c$$$C
 989 c$$$   40 call cartder
 990 c$$$C
 991 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
 992 c$$$C
 993 c$$$
 994 c$$$      ig=0
 995 c$$$      ind=nres-2
 996 c$$$      do i=2,nres-2
 997 c$$$       IF (mask_phi(i+2).eq.1) THEN
 998 c$$$        gphii=0.0D0
 999 c$$$        do j=i+1,nres-1
1000 c$$$          ind=ind+1
1001 c$$$          do k=1,3
1002 c$$$            gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1003 c$$$            gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1004 c$$$          enddo
1005 c$$$        enddo
1006 c$$$        ig=ig+1
1007 c$$$        g(ig)=gphii
1008 c$$$       ELSE
1009 c$$$        ind=ind+nres-1-i
1010 c$$$       ENDIF
1011 c$$$      enddo
1012 c$$$
1013 c$$$
1014 c$$$      ind=0
1015 c$$$      do i=1,nres-2
1016 c$$$       IF (mask_theta(i+2).eq.1) THEN
1017 c$$$        ig=ig+1
1018 c$$$    gthetai=0.0D0
1019 c$$$    do j=i+1,nres-1
1020 c$$$          ind=ind+1
1021 c$$$      do k=1,3
1022 c$$$            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1023 c$$$            gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1024 c$$$          enddo
1025 c$$$        enddo
1026 c$$$        g(ig)=gthetai
1027 c$$$       ELSE
1028 c$$$        ind=ind+nres-1-i
1029 c$$$       ENDIF
1030 c$$$      enddo
1031 c$$$
1032 c$$$      do i=2,nres-1
1033 c$$$    if (itype(i).ne.10) then
1034 c$$$         IF (mask_side(i).eq.1) THEN
1035 c$$$          ig=ig+1
1036 c$$$          galphai=0.0D0
1037 c$$$      do k=1,3
1038 c$$$        galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1039 c$$$          enddo
1040 c$$$          g(ig)=galphai
1041 c$$$         ENDIF
1042 c$$$        endif
1043 c$$$      enddo
1044 c$$$
1045 c$$$
1046 c$$$      do i=2,nres-1
1047 c$$$        if (itype(i).ne.10) then
1048 c$$$         IF (mask_side(i).eq.1) THEN
1049 c$$$          ig=ig+1
1050 c$$$      gomegai=0.0D0
1051 c$$$      do k=1,3
1052 c$$$        gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1053 c$$$          enddo
1054 c$$$      g(ig)=gomegai
1055 c$$$         ENDIF
1056 c$$$        endif
1057 c$$$      enddo
1058 c$$$
1059 c$$$C
1060 c$$$C Add the components corresponding to local energy terms.
1061 c$$$C
1062 c$$$
1063 c$$$      ig=0
1064 c$$$      igall=0
1065 c$$$      do i=4,nres
1066 c$$$        igall=igall+1
1067 c$$$        if (mask_phi(i).eq.1) then
1068 c$$$          ig=ig+1
1069 c$$$          g(ig)=g(ig)+gloc(igall,icg)
1070 c$$$        endif
1071 c$$$      enddo
1072 c$$$
1073 c$$$      do i=3,nres
1074 c$$$        igall=igall+1
1075 c$$$        if (mask_theta(i).eq.1) then
1076 c$$$          ig=ig+1
1077 c$$$          g(ig)=g(ig)+gloc(igall,icg)
1078 c$$$        endif
1079 c$$$      enddo
1080 c$$$
1081 c$$$      do ij=1,2
1082 c$$$      do i=2,nres-1
1083 c$$$        if (itype(i).ne.10) then
1084 c$$$          igall=igall+1
1085 c$$$          if (mask_side(i).eq.1) then
1086 c$$$            ig=ig+1
1087 c$$$            g(ig)=g(ig)+gloc(igall,icg)
1088 c$$$          endif
1089 c$$$        endif
1090 c$$$      enddo
1091 c$$$      enddo
1092 c$$$
1093 c$$$cd      do i=1,ig
1094 c$$$cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1095 c$$$cd      enddo
1096 c$$$
1097 c$$$      return
1098 c$$$      end
1099 c$$$
1100 c$$$C-----------------------------------------------------------------------------
1101 c$$$
1102 c$$$      subroutine etotal_sc(energy_sc)
1103 c$$$      implicit none
1104 c$$$
1105 c$$$c     Includes
1106 c$$$      include 'DIMENSIONS'
1107 c$$$      include 'COMMON.VAR'
1108 c$$$      include 'COMMON.INTERACT'
1109 c$$$      include 'COMMON.DERIV'
1110 c$$$      include 'COMMON.FFIELD'
1111 c$$$
1112 c$$$c     Output arguments
1113 c$$$      double precision energy_sc(0:n_ene)
1114 c$$$
1115 c$$$c     Local variables
1116 c$$$      double precision evdw,escloc
1117 c$$$      integer i,j
1118 c$$$
1119 c$$$
1120 c$$$      do i=1,n_ene
1121 c$$$        energy_sc(i)=0.0D0
1122 c$$$      enddo
1123 c$$$
1124 c$$$      if (mask_r) then
1125 c$$$        call egb_sc(evdw)
1126 c$$$        call esc_sc(escloc)
1127 c$$$      else
1128 c$$$        call egb(evdw)
1129 c$$$        call esc(escloc)
1130 c$$$      endif
1131 c$$$
1132 c$$$      if (evdw.eq.1.0D20) then
1133 c$$$        energy_sc(0)=evdw
1134 c$$$      else
1135 c$$$        energy_sc(0)=wsc*evdw+wscloc*escloc
1136 c$$$      endif
1137 c$$$      energy_sc(1)=evdw
1138 c$$$      energy_sc(12)=escloc
1139 c$$$
1140 c$$$C
1141 c$$$C Sum up the components of the Cartesian gradient.
1142 c$$$C
1143 c$$$      do i=1,nct
1144 c$$$        do j=1,3
1145 c$$$          gradx(j,i,icg)=wsc*gvdwx(j,i)
1146 c$$$        enddo
1147 c$$$      enddo
1148 c$$$
1149 c$$$      return
1150 c$$$      end
1151 c$$$
1152 c$$$C-----------------------------------------------------------------------------
1153 c$$$
1154 c$$$      subroutine egb_sc(evdw)
1155 c$$$C
1156 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1157 c$$$C assuming the Gay-Berne potential of interaction.
1158 c$$$C
1159 c$$$      implicit real*8 (a-h,o-z)
1160 c$$$      include 'DIMENSIONS'
1161 c$$$      include 'COMMON.GEO'
1162 c$$$      include 'COMMON.VAR'
1163 c$$$      include 'COMMON.LOCAL'
1164 c$$$      include 'COMMON.CHAIN'
1165 c$$$      include 'COMMON.DERIV'
1166 c$$$      include 'COMMON.NAMES'
1167 c$$$      include 'COMMON.INTERACT'
1168 c$$$      include 'COMMON.IOUNITS'
1169 c$$$      include 'COMMON.CALC'
1170 c$$$      include 'COMMON.CONTROL'
1171 c$$$      logical lprn
1172 c$$$      evdw=0.0D0
1173 c$$$      energy_dec=.false.
1174 c$$$c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1175 c$$$      evdw=0.0D0
1176 c$$$      lprn=.false.
1177 c$$$c     if (icall.eq.0) lprn=.false.
1178 c$$$      ind=0
1179 c$$$      do i=iatsc_s,iatsc_e
1180 c$$$        itypi=itype(i)
1181 c$$$        itypi1=itype(i+1)
1182 c$$$        xi=c(1,nres+i)
1183 c$$$        yi=c(2,nres+i)
1184 c$$$        zi=c(3,nres+i)
1185 c$$$        dxi=dc_norm(1,nres+i)
1186 c$$$        dyi=dc_norm(2,nres+i)
1187 c$$$        dzi=dc_norm(3,nres+i)
1188 c$$$c        dsci_inv=dsc_inv(itypi)
1189 c$$$        dsci_inv=vbld_inv(i+nres)
1190 c$$$c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1191 c$$$c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1192 c$$$C
1193 c$$$C Calculate SC interaction energy.
1194 c$$$C
1195 c$$$        do iint=1,nint_gr(i)
1196 c$$$          do j=istart(i,iint),iend(i,iint)
1197 c$$$          IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1198 c$$$            ind=ind+1
1199 c$$$            itypj=itype(j)
1200 c$$$c            dscj_inv=dsc_inv(itypj)
1201 c$$$            dscj_inv=vbld_inv(j+nres)
1202 c$$$c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1203 c$$$c     &       1.0d0/vbld(j+nres)
1204 c$$$c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1205 c$$$            sig0ij=sigma(itypi,itypj)
1206 c$$$            chi1=chi(itypi,itypj)
1207 c$$$            chi2=chi(itypj,itypi)
1208 c$$$            chi12=chi1*chi2
1209 c$$$            chip1=chip(itypi)
1210 c$$$            chip2=chip(itypj)
1211 c$$$            chip12=chip1*chip2
1212 c$$$            alf1=alp(itypi)
1213 c$$$            alf2=alp(itypj)
1214 c$$$            alf12=0.5D0*(alf1+alf2)
1215 c$$$C For diagnostics only!!!
1216 c$$$c           chi1=0.0D0
1217 c$$$c           chi2=0.0D0
1218 c$$$c           chi12=0.0D0
1219 c$$$c           chip1=0.0D0
1220 c$$$c           chip2=0.0D0
1221 c$$$c           chip12=0.0D0
1222 c$$$c           alf1=0.0D0
1223 c$$$c           alf2=0.0D0
1224 c$$$c           alf12=0.0D0
1225 c$$$            xj=c(1,nres+j)-xi
1226 c$$$            yj=c(2,nres+j)-yi
1227 c$$$            zj=c(3,nres+j)-zi
1228 c$$$            dxj=dc_norm(1,nres+j)
1229 c$$$            dyj=dc_norm(2,nres+j)
1230 c$$$            dzj=dc_norm(3,nres+j)
1231 c$$$c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1232 c$$$c            write (iout,*) "j",j," dc_norm",
1233 c$$$c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1234 c$$$            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1235 c$$$            rij=dsqrt(rrij)
1236 c$$$C Calculate angle-dependent terms of energy and contributions to their
1237 c$$$C derivatives.
1238 c$$$            call sc_angular
1239 c$$$            sigsq=1.0D0/sigsq
1240 c$$$            sig=sig0ij*dsqrt(sigsq)
1241 c$$$            rij_shift=1.0D0/rij-sig+sig0ij
1242 c$$$c for diagnostics; uncomment
1243 c$$$c            rij_shift=1.2*sig0ij
1244 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1245 c$$$            if (rij_shift.le.0.0D0) then
1246 c$$$              evdw=1.0D20
1247 c$$$cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1248 c$$$cd     &        restyp(itypi),i,restyp(itypj),j,
1249 c$$$cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1250 c$$$              return
1251 c$$$            endif
1252 c$$$            sigder=-sig*sigsq
1253 c$$$c---------------------------------------------------------------
1254 c$$$            rij_shift=1.0D0/rij_shift
1255 c$$$            fac=rij_shift**expon
1256 c$$$            e1=fac*fac*aa(itypi,itypj)
1257 c$$$            e2=fac*bb(itypi,itypj)
1258 c$$$            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1259 c$$$            eps2der=evdwij*eps3rt
1260 c$$$            eps3der=evdwij*eps2rt
1261 c$$$c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1262 c$$$c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1263 c$$$            evdwij=evdwij*eps2rt*eps3rt
1264 c$$$            evdw=evdw+evdwij
1265 c$$$            if (lprn) then
1266 c$$$            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1267 c$$$            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1268 c$$$            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1269 c$$$     &        restyp(itypi),i,restyp(itypj),j,
1270 c$$$     &        epsi,sigm,chi1,chi2,chip1,chip2,
1271 c$$$     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1272 c$$$     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1273 c$$$     &        evdwij
1274 c$$$            endif
1275 c$$$
1276 c$$$            if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1277 c$$$     &                        'evdw',i,j,evdwij
1278 c$$$
1279 c$$$C Calculate gradient components.
1280 c$$$            e1=e1*eps1*eps2rt**2*eps3rt**2
1281 c$$$            fac=-expon*(e1+evdwij)*rij_shift
1282 c$$$            sigder=fac*sigder
1283 c$$$            fac=rij*fac
1284 c$$$c            fac=0.0d0
1285 c$$$C Calculate the radial part of the gradient
1286 c$$$            gg(1)=xj*fac
1287 c$$$            gg(2)=yj*fac
1288 c$$$            gg(3)=zj*fac
1289 c$$$C Calculate angular part of the gradient.
1290 c$$$            call sc_grad
1291 c$$$          ENDIF
1292 c$$$          enddo      ! j
1293 c$$$        enddo        ! iint
1294 c$$$      enddo          ! i
1295 c$$$      energy_dec=.false.
1296 c$$$      return
1297 c$$$      end
1298 c$$$
1299 c$$$c-----------------------------------------------------------------------------
1300 c$$$
1301 c$$$      subroutine esc_sc(escloc)
1302 c$$$C Calculate the local energy of a side chain and its derivatives in the
1303 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1304 c$$$C ALPHA and OMEGA.
1305 c$$$      implicit real*8 (a-h,o-z)
1306 c$$$      include 'DIMENSIONS'
1307 c$$$      include 'COMMON.GEO'
1308 c$$$      include 'COMMON.LOCAL'
1309 c$$$      include 'COMMON.VAR'
1310 c$$$      include 'COMMON.INTERACT'
1311 c$$$      include 'COMMON.DERIV'
1312 c$$$      include 'COMMON.CHAIN'
1313 c$$$      include 'COMMON.IOUNITS'
1314 c$$$      include 'COMMON.NAMES'
1315 c$$$      include 'COMMON.FFIELD'
1316 c$$$      include 'COMMON.CONTROL'
1317 c$$$      double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1318 c$$$     &     ddersc0(3),ddummy(3),xtemp(3),temp(3)
1319 c$$$      common /sccalc/ time11,time12,time112,theti,it,nlobit
1320 c$$$      delta=0.02d0*pi
1321 c$$$      escloc=0.0D0
1322 c$$$c     write (iout,'(a)') 'ESC'
1323 c$$$      do i=loc_start,loc_end
1324 c$$$      IF (mask_side(i).eq.1) THEN
1325 c$$$        it=itype(i)
1326 c$$$        if (it.eq.10) goto 1
1327 c$$$        nlobit=nlob(it)
1328 c$$$c       print *,'i=',i,' it=',it,' nlobit=',nlobit
1329 c$$$c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1330 c$$$        theti=theta(i+1)-pipol
1331 c$$$        x(1)=dtan(theti)
1332 c$$$        x(2)=alph(i)
1333 c$$$        x(3)=omeg(i)
1334 c$$$
1335 c$$$        if (x(2).gt.pi-delta) then
1336 c$$$          xtemp(1)=x(1)
1337 c$$$          xtemp(2)=pi-delta
1338 c$$$          xtemp(3)=x(3)
1339 c$$$          call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1340 c$$$          xtemp(2)=pi
1341 c$$$          call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1342 c$$$          call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1343 c$$$     &        escloci,dersc(2))
1344 c$$$          call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1345 c$$$     &        ddersc0(1),dersc(1))
1346 c$$$          call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1347 c$$$     &        ddersc0(3),dersc(3))
1348 c$$$          xtemp(2)=pi-delta
1349 c$$$          call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1350 c$$$          xtemp(2)=pi
1351 c$$$          call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1352 c$$$          call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1353 c$$$     &            dersc0(2),esclocbi,dersc02)
1354 c$$$          call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1355 c$$$     &            dersc12,dersc01)
1356 c$$$          call splinthet(x(2),0.5d0*delta,ss,ssd)
1357 c$$$          dersc0(1)=dersc01
1358 c$$$          dersc0(2)=dersc02
1359 c$$$          dersc0(3)=0.0d0
1360 c$$$          do k=1,3
1361 c$$$            dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1362 c$$$          enddo
1363 c$$$          dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1364 c$$$c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1365 c$$$c    &             esclocbi,ss,ssd
1366 c$$$          escloci=ss*escloci+(1.0d0-ss)*esclocbi
1367 c$$$c         escloci=esclocbi
1368 c$$$c         write (iout,*) escloci
1369 c$$$        else if (x(2).lt.delta) then
1370 c$$$          xtemp(1)=x(1)
1371 c$$$          xtemp(2)=delta
1372 c$$$          xtemp(3)=x(3)
1373 c$$$          call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1374 c$$$          xtemp(2)=0.0d0
1375 c$$$          call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1376 c$$$          call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1377 c$$$     &        escloci,dersc(2))
1378 c$$$          call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1379 c$$$     &        ddersc0(1),dersc(1))
1380 c$$$          call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1381 c$$$     &        ddersc0(3),dersc(3))
1382 c$$$          xtemp(2)=delta
1383 c$$$          call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1384 c$$$          xtemp(2)=0.0d0
1385 c$$$          call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1386 c$$$          call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1387 c$$$     &            dersc0(2),esclocbi,dersc02)
1388 c$$$          call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1389 c$$$     &            dersc12,dersc01)
1390 c$$$          dersc0(1)=dersc01
1391 c$$$          dersc0(2)=dersc02
1392 c$$$          dersc0(3)=0.0d0
1393 c$$$          call splinthet(x(2),0.5d0*delta,ss,ssd)
1394 c$$$          do k=1,3
1395 c$$$            dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1396 c$$$          enddo
1397 c$$$          dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1398 c$$$c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1399 c$$$c    &             esclocbi,ss,ssd
1400 c$$$          escloci=ss*escloci+(1.0d0-ss)*esclocbi
1401 c$$$c         write (iout,*) escloci
1402 c$$$        else
1403 c$$$          call enesc(x,escloci,dersc,ddummy,.false.)
1404 c$$$        endif
1405 c$$$
1406 c$$$        escloc=escloc+escloci
1407 c$$$        if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1408 c$$$     &     'escloc',i,escloci
1409 c$$$c       write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1410 c$$$
1411 c$$$        gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1412 c$$$     &   wscloc*dersc(1)
1413 c$$$        gloc(ialph(i,1),icg)=wscloc*dersc(2)
1414 c$$$        gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1415 c$$$    1   continue
1416 c$$$      ENDIF
1417 c$$$      enddo
1418 c$$$      return
1419 c$$$      end
1420 c$$$
1421 c$$$C-----------------------------------------------------------------------------
1422 c$$$
1423 c$$$      subroutine egb_ij(i_sc,j_sc,evdw)
1424 c$$$C
1425 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1426 c$$$C assuming the Gay-Berne potential of interaction.
1427 c$$$C
1428 c$$$      implicit real*8 (a-h,o-z)
1429 c$$$      include 'DIMENSIONS'
1430 c$$$      include 'COMMON.GEO'
1431 c$$$      include 'COMMON.VAR'
1432 c$$$      include 'COMMON.LOCAL'
1433 c$$$      include 'COMMON.CHAIN'
1434 c$$$      include 'COMMON.DERIV'
1435 c$$$      include 'COMMON.NAMES'
1436 c$$$      include 'COMMON.INTERACT'
1437 c$$$      include 'COMMON.IOUNITS'
1438 c$$$      include 'COMMON.CALC'
1439 c$$$      include 'COMMON.CONTROL'
1440 c$$$      logical lprn
1441 c$$$      evdw=0.0D0
1442 c$$$      energy_dec=.false.
1443 c$$$c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1444 c$$$      evdw=0.0D0
1445 c$$$      lprn=.false.
1446 c$$$      ind=0
1447 c$$$c$$$      do i=iatsc_s,iatsc_e
1448 c$$$      i=i_sc
1449 c$$$        itypi=itype(i)
1450 c$$$        itypi1=itype(i+1)
1451 c$$$        xi=c(1,nres+i)
1452 c$$$        yi=c(2,nres+i)
1453 c$$$        zi=c(3,nres+i)
1454 c$$$        dxi=dc_norm(1,nres+i)
1455 c$$$        dyi=dc_norm(2,nres+i)
1456 c$$$        dzi=dc_norm(3,nres+i)
1457 c$$$c        dsci_inv=dsc_inv(itypi)
1458 c$$$        dsci_inv=vbld_inv(i+nres)
1459 c$$$c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1460 c$$$c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1461 c$$$C
1462 c$$$C Calculate SC interaction energy.
1463 c$$$C
1464 c$$$c$$$        do iint=1,nint_gr(i)
1465 c$$$c$$$          do j=istart(i,iint),iend(i,iint)
1466 c$$$        j=j_sc
1467 c$$$            ind=ind+1
1468 c$$$            itypj=itype(j)
1469 c$$$c            dscj_inv=dsc_inv(itypj)
1470 c$$$            dscj_inv=vbld_inv(j+nres)
1471 c$$$c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1472 c$$$c     &       1.0d0/vbld(j+nres)
1473 c$$$c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1474 c$$$            sig0ij=sigma(itypi,itypj)
1475 c$$$            chi1=chi(itypi,itypj)
1476 c$$$            chi2=chi(itypj,itypi)
1477 c$$$            chi12=chi1*chi2
1478 c$$$            chip1=chip(itypi)
1479 c$$$            chip2=chip(itypj)
1480 c$$$            chip12=chip1*chip2
1481 c$$$            alf1=alp(itypi)
1482 c$$$            alf2=alp(itypj)
1483 c$$$            alf12=0.5D0*(alf1+alf2)
1484 c$$$C For diagnostics only!!!
1485 c$$$c           chi1=0.0D0
1486 c$$$c           chi2=0.0D0
1487 c$$$c           chi12=0.0D0
1488 c$$$c           chip1=0.0D0
1489 c$$$c           chip2=0.0D0
1490 c$$$c           chip12=0.0D0
1491 c$$$c           alf1=0.0D0
1492 c$$$c           alf2=0.0D0
1493 c$$$c           alf12=0.0D0
1494 c$$$            xj=c(1,nres+j)-xi
1495 c$$$            yj=c(2,nres+j)-yi
1496 c$$$            zj=c(3,nres+j)-zi
1497 c$$$            dxj=dc_norm(1,nres+j)
1498 c$$$            dyj=dc_norm(2,nres+j)
1499 c$$$            dzj=dc_norm(3,nres+j)
1500 c$$$c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1501 c$$$c            write (iout,*) "j",j," dc_norm",
1502 c$$$c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1503 c$$$            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1504 c$$$            rij=dsqrt(rrij)
1505 c$$$C Calculate angle-dependent terms of energy and contributions to their
1506 c$$$C derivatives.
1507 c$$$            call sc_angular
1508 c$$$            sigsq=1.0D0/sigsq
1509 c$$$            sig=sig0ij*dsqrt(sigsq)
1510 c$$$            rij_shift=1.0D0/rij-sig+sig0ij
1511 c$$$c for diagnostics; uncomment
1512 c$$$c            rij_shift=1.2*sig0ij
1513 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1514 c$$$            if (rij_shift.le.0.0D0) then
1515 c$$$              evdw=1.0D20
1516 c$$$cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1517 c$$$cd     &        restyp(itypi),i,restyp(itypj),j,
1518 c$$$cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1519 c$$$              return
1520 c$$$            endif
1521 c$$$            sigder=-sig*sigsq
1522 c$$$c---------------------------------------------------------------
1523 c$$$            rij_shift=1.0D0/rij_shift
1524 c$$$            fac=rij_shift**expon
1525 c$$$            e1=fac*fac*aa(itypi,itypj)
1526 c$$$            e2=fac*bb(itypi,itypj)
1527 c$$$            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1528 c$$$            eps2der=evdwij*eps3rt
1529 c$$$            eps3der=evdwij*eps2rt
1530 c$$$c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1531 c$$$c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1532 c$$$            evdwij=evdwij*eps2rt*eps3rt
1533 c$$$            evdw=evdw+evdwij
1534 c$$$            if (lprn) then
1535 c$$$            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1536 c$$$            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1537 c$$$            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1538 c$$$     &        restyp(itypi),i,restyp(itypj),j,
1539 c$$$     &        epsi,sigm,chi1,chi2,chip1,chip2,
1540 c$$$     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1541 c$$$     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1542 c$$$     &        evdwij
1543 c$$$            endif
1544 c$$$
1545 c$$$            if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1546 c$$$     &                        'evdw',i,j,evdwij
1547 c$$$
1548 c$$$C Calculate gradient components.
1549 c$$$            e1=e1*eps1*eps2rt**2*eps3rt**2
1550 c$$$            fac=-expon*(e1+evdwij)*rij_shift
1551 c$$$            sigder=fac*sigder
1552 c$$$            fac=rij*fac
1553 c$$$c            fac=0.0d0
1554 c$$$C Calculate the radial part of the gradient
1555 c$$$            gg(1)=xj*fac
1556 c$$$            gg(2)=yj*fac
1557 c$$$            gg(3)=zj*fac
1558 c$$$C Calculate angular part of the gradient.
1559 c$$$            call sc_grad
1560 c$$$c$$$          enddo      ! j
1561 c$$$c$$$        enddo        ! iint
1562 c$$$c$$$      enddo          ! i
1563 c$$$      energy_dec=.false.
1564 c$$$      return
1565 c$$$      end
1566 c$$$
1567 c$$$C-----------------------------------------------------------------------------
1568 c$$$
1569 c$$$      subroutine perturb_side_chain(i,angle)
1570 c$$$      implicit none
1571 c$$$
1572 c$$$c     Includes
1573 c$$$      include 'DIMENSIONS'
1574 c$$$      include 'COMMON.CHAIN'
1575 c$$$      include 'COMMON.GEO'
1576 c$$$      include 'COMMON.VAR'
1577 c$$$      include 'COMMON.LOCAL'
1578 c$$$      include 'COMMON.IOUNITS'
1579 c$$$
1580 c$$$c     External functions
1581 c$$$      external ran_number
1582 c$$$      double precision ran_number
1583 c$$$
1584 c$$$c     Input arguments
1585 c$$$      integer i
1586 c$$$      double precision angle    ! In degrees
1587 c$$$
1588 c$$$c     Local variables
1589 c$$$      integer i_sc
1590 c$$$      double precision rad_ang,rand_v(3),length,cost,sint
1591 c$$$
1592 c$$$
1593 c$$$      i_sc=i+nres
1594 c$$$      rad_ang=angle*deg2rad
1595 c$$$
1596 c$$$      length=0.0
1597 c$$$      do while (length.lt.0.01)
1598 c$$$        rand_v(1)=ran_number(0.01D0,1.0D0)
1599 c$$$        rand_v(2)=ran_number(0.01D0,1.0D0)
1600 c$$$        rand_v(3)=ran_number(0.01D0,1.0D0)
1601 c$$$        length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1602 c$$$     +       rand_v(3)*rand_v(3)
1603 c$$$        length=sqrt(length)
1604 c$$$        rand_v(1)=rand_v(1)/length
1605 c$$$        rand_v(2)=rand_v(2)/length
1606 c$$$        rand_v(3)=rand_v(3)/length
1607 c$$$        cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1608 c$$$     +       rand_v(3)*dc_norm(3,i_sc)
1609 c$$$        length=1.0D0-cost*cost
1610 c$$$        if (length.lt.0.0D0) length=0.0D0
1611 c$$$        length=sqrt(length)
1612 c$$$        rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1613 c$$$        rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1614 c$$$        rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1615 c$$$      enddo
1616 c$$$      rand_v(1)=rand_v(1)/length
1617 c$$$      rand_v(2)=rand_v(2)/length
1618 c$$$      rand_v(3)=rand_v(3)/length
1619 c$$$
1620 c$$$      cost=dcos(rad_ang)
1621 c$$$      sint=dsin(rad_ang)
1622 c$$$      dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1623 c$$$      dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1624 c$$$      dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1625 c$$$      dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1626 c$$$      dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1627 c$$$      dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1628 c$$$      c(1,i_sc)=c(1,i)+dc(1,i_sc)
1629 c$$$      c(2,i_sc)=c(2,i)+dc(2,i_sc)
1630 c$$$      c(3,i_sc)=c(3,i)+dc(3,i_sc)
1631 c$$$
1632 c$$$      call chainbuild_cart
1633 c$$$
1634 c$$$      return
1635 c$$$      end
1636 c$$$
1637 c$$$c----------------------------------------------------------------------------
1638 c$$$
1639 c$$$      subroutine ss_relax3(i_in,j_in)
1640 c$$$      implicit none
1641 c$$$
1642 c$$$c     Includes
1643 c$$$      include 'DIMENSIONS'
1644 c$$$      include 'COMMON.VAR'
1645 c$$$      include 'COMMON.CHAIN'
1646 c$$$      include 'COMMON.IOUNITS'
1647 c$$$      include 'COMMON.INTERACT'
1648 c$$$
1649 c$$$c     External functions
1650 c$$$      external ran_number
1651 c$$$      double precision ran_number
1652 c$$$
1653 c$$$c     Input arguments
1654 c$$$      integer i_in,j_in
1655 c$$$
1656 c$$$c     Local variables
1657 c$$$      double precision energy_sc(0:n_ene),etot
1658 c$$$      double precision org_dc(3),org_dc_norm(3),org_c(3)
1659 c$$$      double precision ang_pert,rand_fact,exp_fact,beta
1660 c$$$      integer n,i_pert,i
1661 c$$$      logical notdone
1662 c$$$
1663 c$$$
1664 c$$$      beta=1.0D0
1665 c$$$
1666 c$$$      mask_r=.true.
1667 c$$$      do i=nnt,nct
1668 c$$$        mask_side(i)=0
1669 c$$$      enddo
1670 c$$$      mask_side(i_in)=1
1671 c$$$      mask_side(j_in)=1
1672 c$$$
1673 c$$$      call etotal_sc(energy_sc)
1674 c$$$      etot=energy_sc(0)
1675 c$$$c      write(iout,'(a,3d15.5)')"     SS_MC_START ",energy_sc(0),
1676 c$$$c     +     energy_sc(1),energy_sc(12)
1677 c$$$
1678 c$$$      notdone=.true.
1679 c$$$      n=0
1680 c$$$      do while (notdone)
1681 c$$$        if (mod(n,2).eq.0) then
1682 c$$$          i_pert=i_in
1683 c$$$        else
1684 c$$$          i_pert=j_in
1685 c$$$        endif
1686 c$$$        n=n+1
1687 c$$$
1688 c$$$        do i=1,3
1689 c$$$          org_dc(i)=dc(i,i_pert+nres)
1690 c$$$          org_dc_norm(i)=dc_norm(i,i_pert+nres)
1691 c$$$          org_c(i)=c(i,i_pert+nres)
1692 c$$$        enddo
1693 c$$$        ang_pert=ran_number(0.0D0,3.0D0)
1694 c$$$        call perturb_side_chain(i_pert,ang_pert)
1695 c$$$        call etotal_sc(energy_sc)
1696 c$$$        exp_fact=exp(beta*(etot-energy_sc(0)))
1697 c$$$        rand_fact=ran_number(0.0D0,1.0D0)
1698 c$$$        if (rand_fact.lt.exp_fact) then
1699 c$$$c          write(iout,'(a,3d15.5)')"     SS_MC_ACCEPT ",energy_sc(0),
1700 c$$$c     +     energy_sc(1),energy_sc(12)
1701 c$$$          etot=energy_sc(0)
1702 c$$$        else
1703 c$$$c          write(iout,'(a,3d15.5)')"     SS_MC_REJECT ",energy_sc(0),
1704 c$$$c     +     energy_sc(1),energy_sc(12)
1705 c$$$          do i=1,3
1706 c$$$            dc(i,i_pert+nres)=org_dc(i)
1707 c$$$            dc_norm(i,i_pert+nres)=org_dc_norm(i)
1708 c$$$            c(i,i_pert+nres)=org_c(i)
1709 c$$$          enddo
1710 c$$$        endif
1711 c$$$
1712 c$$$        if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1713 c$$$      enddo
1714 c$$$
1715 c$$$      mask_r=.false.
1716 c$$$
1717 c$$$      return
1718 c$$$      end
1719 c$$$
1720 c$$$c----------------------------------------------------------------------------
1721 c$$$
1722 c$$$      subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1723 c$$$      implicit none
1724 c$$$      include 'DIMENSIONS'
1725 c$$$      integer liv,lv
1726 c$$$      parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1727 c$$$*********************************************************************
1728 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1729 c$$$* the calling subprogram.                                           *
1730 c$$$* when d(i)=1.0, then v(35) is the length of the initial step,      *
1731 c$$$* calculated in the usual pythagorean way.                          *
1732 c$$$* absolute convergence occurs when the function is within v(31) of  *
1733 c$$$* zero. unless you know the minimum value in advance, abs convg     *
1734 c$$$* is probably not useful.                                           *
1735 c$$$* relative convergence is when the model predicts that the function *
1736 c$$$* will decrease by less than v(32)*abs(fun).                        *
1737 c$$$*********************************************************************
1738 c$$$      include 'COMMON.IOUNITS'
1739 c$$$      include 'COMMON.VAR'
1740 c$$$      include 'COMMON.GEO'
1741 c$$$      include 'COMMON.MINIM'
1742 c$$$      include 'COMMON.CHAIN'
1743 c$$$
1744 c$$$      double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1745 c$$$      common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1746 c$$$     +     orig_ss_dist(maxres2,maxres2)
1747 c$$$
1748 c$$$      double precision etot
1749 c$$$      integer iretcode,nfun,i_in,j_in
1750 c$$$
1751 c$$$      external dist
1752 c$$$      double precision dist
1753 c$$$      external ss_func,fdum
1754 c$$$      double precision ss_func,fdum
1755 c$$$
1756 c$$$      integer iv(liv),uiparm(2)
1757 c$$$      double precision v(lv),x(maxres6),d(maxres6),rdum
1758 c$$$      integer i,j,k
1759 c$$$
1760 c$$$
1761 c$$$      call deflt(2,iv,liv,lv,v)
1762 c$$$* 12 means fresh start, dont call deflt
1763 c$$$      iv(1)=12
1764 c$$$* max num of fun calls
1765 c$$$      if (maxfun.eq.0) maxfun=500
1766 c$$$      iv(17)=maxfun
1767 c$$$* max num of iterations
1768 c$$$      if (maxmin.eq.0) maxmin=1000
1769 c$$$      iv(18)=maxmin
1770 c$$$* controls output
1771 c$$$      iv(19)=2
1772 c$$$* selects output unit
1773 c$$$c      iv(21)=iout
1774 c$$$      iv(21)=0
1775 c$$$* 1 means to print out result
1776 c$$$      iv(22)=0
1777 c$$$* 1 means to print out summary stats
1778 c$$$      iv(23)=0
1779 c$$$* 1 means to print initial x and d
1780 c$$$      iv(24)=0
1781 c$$$* min val for v(radfac) default is 0.1
1782 c$$$      v(24)=0.1D0
1783 c$$$* max val for v(radfac) default is 4.0
1784 c$$$      v(25)=2.0D0
1785 c$$$c     v(25)=4.0D0
1786 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1787 c$$$* the sumsl default is 0.1
1788 c$$$      v(26)=0.1D0
1789 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1790 c$$$* the sumsl default is 100*machep
1791 c$$$      v(34)=v(34)/100.0D0
1792 c$$$* absolute convergence
1793 c$$$      if (tolf.eq.0.0D0) tolf=1.0D-4
1794 c$$$      v(31)=tolf
1795 c$$$      v(31)=1.0D-1
1796 c$$$* relative convergence
1797 c$$$      if (rtolf.eq.0.0D0) rtolf=1.0D-4
1798 c$$$      v(32)=rtolf
1799 c$$$      v(32)=1.0D-1
1800 c$$$* controls initial step size
1801 c$$$      v(35)=1.0D-1
1802 c$$$* large vals of d correspond to small components of step
1803 c$$$      do i=1,6*nres
1804 c$$$        d(i)=1.0D0
1805 c$$$      enddo
1806 c$$$
1807 c$$$      do i=0,2*nres
1808 c$$$        do j=1,3
1809 c$$$          orig_ss_dc(j,i)=dc(j,i)
1810 c$$$        enddo
1811 c$$$      enddo
1812 c$$$      call geom_to_var(nvar,orig_ss_var)
1813 c$$$
1814 c$$$      do i=1,nres
1815 c$$$        do j=i,nres
1816 c$$$          orig_ss_dist(j,i)=dist(j,i)
1817 c$$$          orig_ss_dist(j+nres,i)=dist(j+nres,i)
1818 c$$$          orig_ss_dist(j,i+nres)=dist(j,i+nres)
1819 c$$$          orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1820 c$$$        enddo
1821 c$$$      enddo
1822 c$$$
1823 c$$$      k=0
1824 c$$$      do i=1,nres-1
1825 c$$$        do j=1,3
1826 c$$$          k=k+1
1827 c$$$          x(k)=dc(j,i)
1828 c$$$        enddo
1829 c$$$      enddo
1830 c$$$      do i=2,nres-1
1831 c$$$        if (ialph(i,1).gt.0) then
1832 c$$$        do j=1,3
1833 c$$$          k=k+1
1834 c$$$          x(k)=dc(j,i+nres)
1835 c$$$        enddo
1836 c$$$        endif
1837 c$$$      enddo
1838 c$$$
1839 c$$$      uiparm(1)=i_in
1840 c$$$      uiparm(2)=j_in
1841 c$$$      call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1842 c$$$      etot=v(10)
1843 c$$$      iretcode=iv(1)
1844 c$$$      nfun=iv(6)+iv(30)
1845 c$$$
1846 c$$$      k=0
1847 c$$$      do i=1,nres-1
1848 c$$$        do j=1,3
1849 c$$$          k=k+1
1850 c$$$          dc(j,i)=x(k)
1851 c$$$        enddo
1852 c$$$      enddo
1853 c$$$      do i=2,nres-1
1854 c$$$        if (ialph(i,1).gt.0) then
1855 c$$$        do j=1,3
1856 c$$$          k=k+1
1857 c$$$          dc(j,i+nres)=x(k)
1858 c$$$        enddo
1859 c$$$        endif
1860 c$$$      enddo
1861 c$$$      call chainbuild_cart
1862 c$$$
1863 c$$$      return
1864 c$$$      end
1865 c$$$
1866 c$$$C-----------------------------------------------------------------------------
1867 c$$$
1868 c$$$      subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1869 c$$$      implicit none
1870 c$$$      include 'DIMENSIONS'
1871 c$$$      include 'COMMON.DERIV'
1872 c$$$      include 'COMMON.IOUNITS'
1873 c$$$      include 'COMMON.VAR'
1874 c$$$      include 'COMMON.CHAIN'
1875 c$$$      include 'COMMON.INTERACT'
1876 c$$$      include 'COMMON.SBRIDGE'
1877 c$$$
1878 c$$$      double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1879 c$$$      common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1880 c$$$     +     orig_ss_dist(maxres2,maxres2)
1881 c$$$
1882 c$$$      integer n
1883 c$$$      double precision x(maxres6)
1884 c$$$      integer nf
1885 c$$$      double precision f
1886 c$$$      integer uiparm(2)
1887 c$$$      real*8 urparm(1)
1888 c$$$      external ufparm
1889 c$$$      double precision ufparm
1890 c$$$
1891 c$$$      external dist
1892 c$$$      double precision dist
1893 c$$$
1894 c$$$      integer i,j,k,ss_i,ss_j
1895 c$$$      double precision tempf,var(maxvar)
1896 c$$$
1897 c$$$
1898 c$$$      ss_i=uiparm(1)
1899 c$$$      ss_j=uiparm(2)
1900 c$$$      f=0.0D0
1901 c$$$
1902 c$$$      k=0
1903 c$$$      do i=1,nres-1
1904 c$$$        do j=1,3
1905 c$$$          k=k+1
1906 c$$$          dc(j,i)=x(k)
1907 c$$$        enddo
1908 c$$$      enddo
1909 c$$$      do i=2,nres-1
1910 c$$$        if (ialph(i,1).gt.0) then
1911 c$$$        do j=1,3
1912 c$$$          k=k+1
1913 c$$$          dc(j,i+nres)=x(k)
1914 c$$$        enddo
1915 c$$$        endif
1916 c$$$      enddo
1917 c$$$      call chainbuild_cart
1918 c$$$
1919 c$$$      call geom_to_var(nvar,var)
1920 c$$$
1921 c$$$c     Constraints on all angles
1922 c$$$      do i=1,nvar
1923 c$$$        tempf=var(i)-orig_ss_var(i)
1924 c$$$        f=f+tempf*tempf
1925 c$$$      enddo
1926 c$$$
1927 c$$$c     Constraints on all distances
1928 c$$$      do i=1,nres-1
1929 c$$$        if (i.gt.1) then
1930 c$$$          tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1931 c$$$          f=f+tempf*tempf
1932 c$$$        endif
1933 c$$$        do j=i+1,nres
1934 c$$$          tempf=dist(j,i)-orig_ss_dist(j,i)
1935 c$$$          if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1936 c$$$          tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1937 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
1938 c$$$          tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1939 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
1940 c$$$          tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1941 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
1942 c$$$        enddo
1943 c$$$      enddo
1944 c$$$
1945 c$$$c     Constraints for the relevant CYS-CYS
1946 c$$$      tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1947 c$$$      f=f+tempf*tempf
1948 c$$$CCCCCCCCCCCCCCCCC      ADD SOME ANGULAR STUFF
1949 c$$$
1950 c$$$c$$$      if (nf.ne.nfl) then
1951 c$$$c$$$        write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1952 c$$$c$$$     +       f,dist(5+nres,14+nres)
1953 c$$$c$$$      endif
1954 c$$$
1955 c$$$      nfl=nf
1956 c$$$
1957 c$$$      return
1958 c$$$      end
1959 c$$$
1960 c$$$C-----------------------------------------------------------------------------