source/cluster/wham/src-M/ssMD.F

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