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