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