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