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 #ifdef WHAM
 766       subroutine read_ssHist
 767       implicit none
 768
 769 c     Includes
 770       include 'DIMENSIONS'
 771       include "DIMENSIONS.FREE"
 772       include 'COMMON.FREE'
 773
 774 c     Local variables
 775       integer i,j
 776       character*80 controlcard
 777
 778       do i=1,dyn_nssHist
 779         call card_concat(controlcard,.true.)
 780         read(controlcard,*)
 781      &       dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
 782       enddo
 783
 784       return
 785       end
 786 #endif
 787
 788 c----------------------------------------------------------------------------
 789
 790
 791 C-----------------------------------------------------------------------------
 792 C-----------------------------------------------------------------------------
 793 C-----------------------------------------------------------------------------
 794 C-----------------------------------------------------------------------------
 795 C-----------------------------------------------------------------------------
 796 C-----------------------------------------------------------------------------
 797 C-----------------------------------------------------------------------------
 798
 799 c$$$c-----------------------------------------------------------------------------
 800 c$$$
 801 c$$$      subroutine ss_relax(i_in,j_in)
 802 c$$$      implicit none
 803 c$$$
 804 c$$$c     Includes
 805 c$$$      include 'DIMENSIONS'
 806 c$$$      include 'COMMON.VAR'
 807 c$$$      include 'COMMON.CHAIN'
 808 c$$$      include 'COMMON.IOUNITS'
 809 c$$$      include 'COMMON.INTERACT'
 810 c$$$
 811 c$$$c     Input arguments
 812 c$$$      integer i_in,j_in
 813 c$$$
 814 c$$$c     Local variables
 815 c$$$      integer i,iretcode,nfun_sc
 816 c$$$      logical scfail
 817 c$$$      double precision var(maxvar),e_sc,etot
 818 c$$$
 819 c$$$
 820 c$$$      mask_r=.true.
 821 c$$$      do i=nnt,nct
 822 c$$$        mask_side(i)=0
 823 c$$$      enddo
 824 c$$$      mask_side(i_in)=1
 825 c$$$      mask_side(j_in)=1
 826 c$$$
 827 c$$$c     Minimize the two selected side-chains
 828 c$$$      call overlap_sc(scfail)  ! Better not fail!
 829 c$$$      call minimize_sc(e_sc,var,iretcode,nfun_sc)
 830 c$$$
 831 c$$$      mask_r=.false.
 832 c$$$
 833 c$$$      return
 834 c$$$      end
 835 c$$$
 836 c$$$c-------------------------------------------------------------
 837 c$$$
 838 c$$$      subroutine minimize_sc(etot_sc,iretcode,nfun)
 839 c$$$c     Minimize side-chains only, starting from geom but without modifying
 840 c$$$c     bond lengths.
 841 c$$$c     If mask_r is already set, only the selected side-chains are minimized,
 842 c$$$c     otherwise all side-chains are minimized keeping the backbone frozen.
 843 c$$$      implicit none
 844 c$$$
 845 c$$$c     Includes
 846 c$$$      include 'DIMENSIONS'
 847 c$$$      include 'COMMON.IOUNITS'
 848 c$$$      include 'COMMON.VAR'
 849 c$$$      include 'COMMON.CHAIN'
 850 c$$$      include 'COMMON.GEO'
 851 c$$$      include 'COMMON.MINIM'
 852 c$$$      integer icall
 853 c$$$      common /srutu/ icall
 854 c$$$
 855 c$$$c     Output arguments
 856 c$$$      double precision etot_sc
 857 c$$$      integer iretcode,nfun
 858 c$$$
 859 c$$$c     External functions/subroutines
 860 c$$$      external func_sc,grad_sc,fdum
 861 c$$$
 862 c$$$c     Local variables
 863 c$$$      integer liv,lv
 864 c$$$      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
 865 c$$$      integer iv(liv)
 866 c$$$      double precision rdum(1)
 867 c$$$      double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
 868 c$$$      integer idum(1)
 869 c$$$      integer i,nvar_restr
 870 c$$$
 871 c$$$
 872 c$$$cmc      start_minim=.true.
 873 c$$$      call deflt(2,iv,liv,lv,v)
 874 c$$$* 12 means fresh start, dont call deflt
 875 c$$$      iv(1)=12
 876 c$$$* max num of fun calls
 877 c$$$      if (maxfun.eq.0) maxfun=500
 878 c$$$      iv(17)=maxfun
 879 c$$$* max num of iterations
 880 c$$$      if (maxmin.eq.0) maxmin=1000
 881 c$$$      iv(18)=maxmin
 882 c$$$* controls output
 883 c$$$      iv(19)=1
 884 c$$$* selects output unit
 885 c$$$      iv(21)=0
 886 c$$$c      iv(21)=iout               ! DEBUG
 887 c$$$c      iv(21)=8                  ! DEBUG
 888 c$$$* 1 means to print out result
 889 c$$$      iv(22)=0
 890 c$$$c      iv(22)=1                  ! DEBUG
 891 c$$$* 1 means to print out summary stats
 892 c$$$      iv(23)=0
 893 c$$$c      iv(23)=1                  ! DEBUG
 894 c$$$* 1 means to print initial x and d
 895 c$$$      iv(24)=0
 896 c$$$c      iv(24)=1                  ! DEBUG
 897 c$$$* min val for v(radfac) default is 0.1
 898 c$$$      v(24)=0.1D0
 899 c$$$* max val for v(radfac) default is 4.0
 900 c$$$      v(25)=2.0D0
 901 c$$$c     v(25)=4.0D0
 902 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
 903 c$$$* the sumsl default is 0.1
 904 c$$$      v(26)=0.1D0
 905 c$$$* false conv if (act fnctn decrease) .lt. v(34)
 906 c$$$* the sumsl default is 100*machep
 907 c$$$      v(34)=v(34)/100.0D0
 908 c$$$* absolute convergence
 909 c$$$      if (tolf.eq.0.0D0) tolf=1.0D-4
 910 c$$$      v(31)=tolf
 911 c$$$* relative convergence
 912 c$$$      if (rtolf.eq.0.0D0) rtolf=1.0D-1
 913 c$$$      v(32)=rtolf
 914 c$$$* controls initial step size
 915 c$$$       v(35)=1.0D-1
 916 c$$$* large vals of d correspond to small components of step
 917 c$$$      do i=1,nphi
 918 c$$$        d(i)=1.0D-1
 919 c$$$      enddo
 920 c$$$      do i=nphi+1,nvar
 921 c$$$        d(i)=1.0D-1
 922 c$$$      enddo
 923 c$$$
 924 c$$$      call geom_to_var(nvar,x)
 925 c$$$      IF (mask_r) THEN
 926 c$$$        do i=1,nres             ! Just in case...
 927 c$$$          mask_phi(i)=0
 928 c$$$          mask_theta(i)=0
 929 c$$$        enddo
 930 c$$$        call x2xx(x,xx,nvar_restr)
 931 c$$$        call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
 932 c$$$     &       iv,liv,lv,v,idum,rdum,fdum)
 933 c$$$        call xx2x(x,xx)
 934 c$$$      ELSE
 935 c$$$c     When minimizing ALL side-chains, etotal_sc is a little
 936 c$$$c     faster if we don't set mask_r
 937 c$$$        do i=1,nres
 938 c$$$          mask_phi(i)=0
 939 c$$$          mask_theta(i)=0
 940 c$$$          mask_side(i)=1
 941 c$$$        enddo
 942 c$$$        call x2xx(x,xx,nvar_restr)
 943 c$$$        call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
 944 c$$$     &       iv,liv,lv,v,idum,rdum,fdum)
 945 c$$$        call xx2x(x,xx)
 946 c$$$      ENDIF
 947 c$$$      call var_to_geom(nvar,x)
 948 c$$$      call chainbuild_sc
 949 c$$$      etot_sc=v(10)
 950 c$$$      iretcode=iv(1)
 951 c$$$      nfun=iv(6)
 952 c$$$      return
 953 c$$$      end
 954 c$$$
 955 c$$$C--------------------------------------------------------------------------
 956 c$$$
 957 c$$$      subroutine chainbuild_sc
 958 c$$$      implicit none
 959 c$$$      include 'DIMENSIONS'
 960 c$$$      include 'COMMON.VAR'
 961 c$$$      include 'COMMON.INTERACT'
 962 c$$$
 963 c$$$c     Local variables
 964 c$$$      integer i
 965 c$$$
 966 c$$$
 967 c$$$      do i=nnt,nct
 968 c$$$        if (.not.mask_r .or. mask_side(i).eq.1) then
 969 c$$$          call locate_side_chain(i)
 970 c$$$        endif
 971 c$$$      enddo
 972 c$$$
 973 c$$$      return
 974 c$$$      end
 975 c$$$
 976 c$$$C--------------------------------------------------------------------------
 977 c$$$
 978 c$$$      subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
 979 c$$$      implicit none
 980 c$$$
 981 c$$$c     Includes
 982 c$$$      include 'DIMENSIONS'
 983 c$$$      include 'COMMON.DERIV'
 984 c$$$      include 'COMMON.VAR'
 985 c$$$      include 'COMMON.MINIM'
 986 c$$$      include 'COMMON.IOUNITS'
 987 c$$$
 988 c$$$c     Input arguments
 989 c$$$      integer n
 990 c$$$      double precision x(maxvar)
 991 c$$$      double precision ufparm
 992 c$$$      external ufparm
 993 c$$$
 994 c$$$c     Input/Output arguments
 995 c$$$      integer nf
 996 c$$$      integer uiparm(1)
 997 c$$$      double precision urparm(1)
 998 c$$$
 999 c$$$c     Output arguments
1000 c$$$      double precision f
1001 c$$$
1002 c$$$c     Local variables
1003 c$$$      double precision energia(0:n_ene)
1004 c$$$#ifdef OSF
1005 c$$$c     Variables used to intercept NaNs
1006 c$$$      double precision x_sum
1007 c$$$      integer i_NAN
1008 c$$$#endif
1009 c$$$
1010 c$$$
1011 c$$$      nfl=nf
1012 c$$$      icg=mod(nf,2)+1
1013 c$$$
1014 c$$$#ifdef OSF
1015 c$$$c     Intercept NaNs in the coordinates, before calling etotal_sc
1016 c$$$      x_sum=0.D0
1017 c$$$      do i_NAN=1,n
1018 c$$$        x_sum=x_sum+x(i_NAN)
1019 c$$$      enddo
1020 c$$$c     Calculate the energy only if the coordinates are ok
1021 c$$$      if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
1022 c$$$        write(iout,*)"   *** func_restr_sc : Found NaN in coordinates"
1023 c$$$        f=1.0D+77
1024 c$$$        nf=0
1025 c$$$      else
1026 c$$$#endif
1027 c$$$
1028 c$$$      call var_to_geom_restr(n,x)
1029 c$$$      call zerograd
1030 c$$$      call chainbuild_sc
1031 c$$$      call etotal_sc(energia(0))
1032 c$$$      f=energia(0)
1033 c$$$      if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
1034 c$$$
1035 c$$$#ifdef OSF
1036 c$$$      endif
1037 c$$$#endif
1038 c$$$
1039 c$$$      return
1040 c$$$      end
1041 c$$$
1042 c$$$c-------------------------------------------------------
1043 c$$$
1044 c$$$      subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
1045 c$$$      implicit none
1046 c$$$
1047 c$$$c     Includes
1048 c$$$      include 'DIMENSIONS'
1049 c$$$      include 'COMMON.CHAIN'
1050 c$$$      include 'COMMON.DERIV'
1051 c$$$      include 'COMMON.VAR'
1052 c$$$      include 'COMMON.INTERACT'
1053 c$$$      include 'COMMON.MINIM'
1054 c$$$
1055 c$$$c     Input arguments
1056 c$$$      integer n
1057 c$$$      double precision x(maxvar)
1058 c$$$      double precision ufparm
1059 c$$$      external ufparm
1060 c$$$
1061 c$$$c     Input/Output arguments
1062 c$$$      integer nf
1063 c$$$      integer uiparm(1)
1064 c$$$      double precision urparm(1)
1065 c$$$
1066 c$$$c     Output arguments
1067 c$$$      double precision g(maxvar)
1068 c$$$
1069 c$$$c     Local variables
1070 c$$$      double precision f,gphii,gthetai,galphai,gomegai
1071 c$$$      integer ig,ind,i,j,k,igall,ij
1072 c$$$
1073 c$$$
1074 c$$$      icg=mod(nf,2)+1
1075 c$$$      if (nf-nfl+1) 20,30,40
1076 c$$$   20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1077 c$$$c     write (iout,*) 'grad 20'
1078 c$$$      if (nf.eq.0) return
1079 c$$$      goto 40
1080 c$$$   30 call var_to_geom_restr(n,x)
1081 c$$$      call chainbuild_sc
1082 c$$$C
1083 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1084 c$$$C
1085 c$$$   40 call cartder
1086 c$$$C
1087 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1088 c$$$C
1089 c$$$
1090 c$$$      ig=0
1091 c$$$      ind=nres-2
1092 c$$$      do i=2,nres-2
1093 c$$$       IF (mask_phi(i+2).eq.1) THEN
1094 c$$$        gphii=0.0D0
1095 c$$$        do j=i+1,nres-1
1096 c$$$          ind=ind+1
1097 c$$$          do k=1,3
1098 c$$$            gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1099 c$$$            gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1100 c$$$          enddo
1101 c$$$        enddo
1102 c$$$        ig=ig+1
1103 c$$$        g(ig)=gphii
1104 c$$$       ELSE
1105 c$$$        ind=ind+nres-1-i
1106 c$$$       ENDIF
1107 c$$$      enddo
1108 c$$$
1109 c$$$
1110 c$$$      ind=0
1111 c$$$      do i=1,nres-2
1112 c$$$       IF (mask_theta(i+2).eq.1) THEN
1113 c$$$        ig=ig+1
1114 c$$$    gthetai=0.0D0
1115 c$$$    do j=i+1,nres-1
1116 c$$$          ind=ind+1
1117 c$$$      do k=1,3
1118 c$$$            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1119 c$$$            gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1120 c$$$          enddo
1121 c$$$        enddo
1122 c$$$        g(ig)=gthetai
1123 c$$$       ELSE
1124 c$$$        ind=ind+nres-1-i
1125 c$$$       ENDIF
1126 c$$$      enddo
1127 c$$$
1128 c$$$      do i=2,nres-1
1129 c$$$    if (itype(i).ne.10) then
1130 c$$$         IF (mask_side(i).eq.1) THEN
1131 c$$$          ig=ig+1
1132 c$$$          galphai=0.0D0
1133 c$$$      do k=1,3
1134 c$$$        galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1135 c$$$          enddo
1136 c$$$          g(ig)=galphai
1137 c$$$         ENDIF
1138 c$$$        endif
1139 c$$$      enddo
1140 c$$$
1141 c$$$
1142 c$$$      do i=2,nres-1
1143 c$$$        if (itype(i).ne.10) then
1144 c$$$         IF (mask_side(i).eq.1) THEN
1145 c$$$          ig=ig+1
1146 c$$$      gomegai=0.0D0
1147 c$$$      do k=1,3
1148 c$$$        gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1149 c$$$          enddo
1150 c$$$      g(ig)=gomegai
1151 c$$$         ENDIF
1152 c$$$        endif
1153 c$$$      enddo
1154 c$$$
1155 c$$$C
1156 c$$$C Add the components corresponding to local energy terms.
1157 c$$$C
1158 c$$$
1159 c$$$      ig=0
1160 c$$$      igall=0
1161 c$$$      do i=4,nres
1162 c$$$        igall=igall+1
1163 c$$$        if (mask_phi(i).eq.1) then
1164 c$$$          ig=ig+1
1165 c$$$          g(ig)=g(ig)+gloc(igall,icg)
1166 c$$$        endif
1167 c$$$      enddo
1168 c$$$
1169 c$$$      do i=3,nres
1170 c$$$        igall=igall+1
1171 c$$$        if (mask_theta(i).eq.1) then
1172 c$$$          ig=ig+1
1173 c$$$          g(ig)=g(ig)+gloc(igall,icg)
1174 c$$$        endif
1175 c$$$      enddo
1176 c$$$
1177 c$$$      do ij=1,2
1178 c$$$      do i=2,nres-1
1179 c$$$        if (itype(i).ne.10) then
1180 c$$$          igall=igall+1
1181 c$$$          if (mask_side(i).eq.1) then
1182 c$$$            ig=ig+1
1183 c$$$            g(ig)=g(ig)+gloc(igall,icg)
1184 c$$$          endif
1185 c$$$        endif
1186 c$$$      enddo
1187 c$$$      enddo
1188 c$$$
1189 c$$$cd      do i=1,ig
1190 c$$$cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1191 c$$$cd      enddo
1192 c$$$
1193 c$$$      return
1194 c$$$      end
1195 c$$$
1196 c$$$C-----------------------------------------------------------------------------
1197 c$$$
1198 c$$$      subroutine etotal_sc(energy_sc)
1199 c$$$      implicit none
1200 c$$$
1201 c$$$c     Includes
1202 c$$$      include 'DIMENSIONS'
1203 c$$$      include 'COMMON.VAR'
1204 c$$$      include 'COMMON.INTERACT'
1205 c$$$      include 'COMMON.DERIV'
1206 c$$$      include 'COMMON.FFIELD'
1207 c$$$
1208 c$$$c     Output arguments
1209 c$$$      double precision energy_sc(0:n_ene)
1210 c$$$
1211 c$$$c     Local variables
1212 c$$$      double precision evdw,escloc
1213 c$$$      integer i,j
1214 c$$$
1215 c$$$
1216 c$$$      do i=1,n_ene
1217 c$$$        energy_sc(i)=0.0D0
1218 c$$$      enddo
1219 c$$$
1220 c$$$      if (mask_r) then
1221 c$$$        call egb_sc(evdw)
1222 c$$$        call esc_sc(escloc)
1223 c$$$      else
1224 c$$$        call egb(evdw)
1225 c$$$        call esc(escloc)
1226 c$$$      endif
1227 c$$$
1228 c$$$      if (evdw.eq.1.0D20) then
1229 c$$$        energy_sc(0)=evdw
1230 c$$$      else
1231 c$$$        energy_sc(0)=wsc*evdw+wscloc*escloc
1232 c$$$      endif
1233 c$$$      energy_sc(1)=evdw
1234 c$$$      energy_sc(12)=escloc
1235 c$$$
1236 c$$$C
1237 c$$$C Sum up the components of the Cartesian gradient.
1238 c$$$C
1239 c$$$      do i=1,nct
1240 c$$$        do j=1,3
1241 c$$$          gradx(j,i,icg)=wsc*gvdwx(j,i)
1242 c$$$        enddo
1243 c$$$      enddo
1244 c$$$
1245 c$$$      return
1246 c$$$      end
1247 c$$$
1248 c$$$C-----------------------------------------------------------------------------
1249 c$$$
1250 c$$$      subroutine egb_sc(evdw)
1251 c$$$C
1252 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1253 c$$$C assuming the Gay-Berne potential of interaction.
1254 c$$$C
1255 c$$$      implicit real*8 (a-h,o-z)
1256 c$$$      include 'DIMENSIONS'
1257 c$$$      include 'COMMON.GEO'
1258 c$$$      include 'COMMON.VAR'
1259 c$$$      include 'COMMON.LOCAL'
1260 c$$$      include 'COMMON.CHAIN'
1261 c$$$      include 'COMMON.DERIV'
1262 c$$$      include 'COMMON.NAMES'
1263 c$$$      include 'COMMON.INTERACT'
1264 c$$$      include 'COMMON.IOUNITS'
1265 c$$$      include 'COMMON.CALC'
1266 c$$$      include 'COMMON.CONTROL'
1267 c$$$      logical lprn
1268 c$$$      evdw=0.0D0
1269 c$$$      energy_dec=.false.
1270 c$$$c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1271 c$$$      evdw=0.0D0
1272 c$$$      lprn=.false.
1273 c$$$c     if (icall.eq.0) lprn=.false.
1274 c$$$      ind=0
1275 c$$$      do i=iatsc_s,iatsc_e
1276 c$$$        itypi=itype(i)
1277 c$$$        itypi1=itype(i+1)
1278 c$$$        xi=c(1,nres+i)
1279 c$$$        yi=c(2,nres+i)
1280 c$$$        zi=c(3,nres+i)
1281 c$$$        dxi=dc_norm(1,nres+i)
1282 c$$$        dyi=dc_norm(2,nres+i)
1283 c$$$        dzi=dc_norm(3,nres+i)
1284 c$$$c        dsci_inv=dsc_inv(itypi)
1285 c$$$        dsci_inv=vbld_inv(i+nres)
1286 c$$$c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1287 c$$$c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1288 c$$$C
1289 c$$$C Calculate SC interaction energy.
1290 c$$$C
1291 c$$$        do iint=1,nint_gr(i)
1292 c$$$          do j=istart(i,iint),iend(i,iint)
1293 c$$$          IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1294 c$$$            ind=ind+1
1295 c$$$            itypj=itype(j)
1296 c$$$c            dscj_inv=dsc_inv(itypj)
1297 c$$$            dscj_inv=vbld_inv(j+nres)
1298 c$$$c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1299 c$$$c     &       1.0d0/vbld(j+nres)
1300 c$$$c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1301 c$$$            sig0ij=sigma(itypi,itypj)
1302 c$$$            chi1=chi(itypi,itypj)
1303 c$$$            chi2=chi(itypj,itypi)
1304 c$$$            chi12=chi1*chi2
1305 c$$$            chip1=chip(itypi)
1306 c$$$            chip2=chip(itypj)
1307 c$$$            chip12=chip1*chip2
1308 c$$$            alf1=alp(itypi)
1309 c$$$            alf2=alp(itypj)
1310 c$$$            alf12=0.5D0*(alf1+alf2)
1311 c$$$C For diagnostics only!!!
1312 c$$$c           chi1=0.0D0
1313 c$$$c           chi2=0.0D0
1314 c$$$c           chi12=0.0D0
1315 c$$$c           chip1=0.0D0
1316 c$$$c           chip2=0.0D0
1317 c$$$c           chip12=0.0D0
1318 c$$$c           alf1=0.0D0
1319 c$$$c           alf2=0.0D0
1320 c$$$c           alf12=0.0D0
1321 c$$$            xj=c(1,nres+j)-xi
1322 c$$$            yj=c(2,nres+j)-yi
1323 c$$$            zj=c(3,nres+j)-zi
1324 c$$$            dxj=dc_norm(1,nres+j)
1325 c$$$            dyj=dc_norm(2,nres+j)
1326 c$$$            dzj=dc_norm(3,nres+j)
1327 c$$$c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1328 c$$$c            write (iout,*) "j",j," dc_norm",
1329 c$$$c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1330 c$$$            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1331 c$$$            rij=dsqrt(rrij)
1332 c$$$C Calculate angle-dependent terms of energy and contributions to their
1333 c$$$C derivatives.
1334 c$$$            call sc_angular
1335 c$$$            sigsq=1.0D0/sigsq
1336 c$$$            sig=sig0ij*dsqrt(sigsq)
1337 c$$$            rij_shift=1.0D0/rij-sig+sig0ij
1338 c$$$c for diagnostics; uncomment
1339 c$$$c            rij_shift=1.2*sig0ij
1340 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1341 c$$$            if (rij_shift.le.0.0D0) then
1342 c$$$              evdw=1.0D20
1343 c$$$cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1344 c$$$cd     &        restyp(itypi),i,restyp(itypj),j,
1345 c$$$cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1346 c$$$              return
1347 c$$$            endif
1348 c$$$            sigder=-sig*sigsq
1349 c$$$c---------------------------------------------------------------
1350 c$$$            rij_shift=1.0D0/rij_shift
1351 c$$$            fac=rij_shift**expon
1352 c$$$            e1=fac*fac*aa(itypi,itypj)
1353 c$$$            e2=fac*bb(itypi,itypj)
1354 c$$$            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1355 c$$$            eps2der=evdwij*eps3rt
1356 c$$$            eps3der=evdwij*eps2rt
1357 c$$$c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1358 c$$$c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1359 c$$$            evdwij=evdwij*eps2rt*eps3rt
1360 c$$$            evdw=evdw+evdwij
1361 c$$$            if (lprn) then
1362 c$$$            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1363 c$$$            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1364 c$$$            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1365 c$$$     &        restyp(itypi),i,restyp(itypj),j,
1366 c$$$     &        epsi,sigm,chi1,chi2,chip1,chip2,
1367 c$$$     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1368 c$$$     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1369 c$$$     &        evdwij
1370 c$$$            endif
1371 c$$$
1372 c$$$            if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1373 c$$$     &                        'evdw',i,j,evdwij
1374 c$$$
1375 c$$$C Calculate gradient components.
1376 c$$$            e1=e1*eps1*eps2rt**2*eps3rt**2
1377 c$$$            fac=-expon*(e1+evdwij)*rij_shift
1378 c$$$            sigder=fac*sigder
1379 c$$$            fac=rij*fac
1380 c$$$c            fac=0.0d0
1381 c$$$C Calculate the radial part of the gradient
1382 c$$$            gg(1)=xj*fac
1383 c$$$            gg(2)=yj*fac
1384 c$$$            gg(3)=zj*fac
1385 c$$$C Calculate angular part of the gradient.
1386 c$$$            call sc_grad
1387 c$$$          ENDIF
1388 c$$$          enddo      ! j
1389 c$$$        enddo        ! iint
1390 c$$$      enddo          ! i
1391 c$$$      energy_dec=.false.
1392 c$$$      return
1393 c$$$      end
1394 c$$$
1395 c$$$c-----------------------------------------------------------------------------
1396 c$$$
1397 c$$$      subroutine esc_sc(escloc)
1398 c$$$C Calculate the local energy of a side chain and its derivatives in the
1399 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1400 c$$$C ALPHA and OMEGA.
1401 c$$$      implicit real*8 (a-h,o-z)
1402 c$$$      include 'DIMENSIONS'
1403 c$$$      include 'COMMON.GEO'
1404 c$$$      include 'COMMON.LOCAL'
1405 c$$$      include 'COMMON.VAR'
1406 c$$$      include 'COMMON.INTERACT'
1407 c$$$      include 'COMMON.DERIV'
1408 c$$$      include 'COMMON.CHAIN'
1409 c$$$      include 'COMMON.IOUNITS'
1410 c$$$      include 'COMMON.NAMES'
1411 c$$$      include 'COMMON.FFIELD'
1412 c$$$      include 'COMMON.CONTROL'
1413 c$$$      double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1414 c$$$     &     ddersc0(3),ddummy(3),xtemp(3),temp(3)
1415 c$$$      common /sccalc/ time11,time12,time112,theti,it,nlobit
1416 c$$$      delta=0.02d0*pi
1417 c$$$      escloc=0.0D0
1418 c$$$c     write (iout,'(a)') 'ESC'
1419 c$$$      do i=loc_start,loc_end
1420 c$$$      IF (mask_side(i).eq.1) THEN
1421 c$$$        it=itype(i)
1422 c$$$        if (it.eq.10) goto 1
1423 c$$$        nlobit=nlob(it)
1424 c$$$c       print *,'i=',i,' it=',it,' nlobit=',nlobit
1425 c$$$c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1426 c$$$        theti=theta(i+1)-pipol
1427 c$$$        x(1)=dtan(theti)
1428 c$$$        x(2)=alph(i)
1429 c$$$        x(3)=omeg(i)
1430 c$$$
1431 c$$$        if (x(2).gt.pi-delta) then
1432 c$$$          xtemp(1)=x(1)
1433 c$$$          xtemp(2)=pi-delta
1434 c$$$          xtemp(3)=x(3)
1435 c$$$          call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1436 c$$$          xtemp(2)=pi
1437 c$$$          call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1438 c$$$          call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1439 c$$$     &        escloci,dersc(2))
1440 c$$$          call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1441 c$$$     &        ddersc0(1),dersc(1))
1442 c$$$          call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1443 c$$$     &        ddersc0(3),dersc(3))
1444 c$$$          xtemp(2)=pi-delta
1445 c$$$          call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1446 c$$$          xtemp(2)=pi
1447 c$$$          call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1448 c$$$          call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1449 c$$$     &            dersc0(2),esclocbi,dersc02)
1450 c$$$          call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1451 c$$$     &            dersc12,dersc01)
1452 c$$$          call splinthet(x(2),0.5d0*delta,ss,ssd)
1453 c$$$          dersc0(1)=dersc01
1454 c$$$          dersc0(2)=dersc02
1455 c$$$          dersc0(3)=0.0d0
1456 c$$$          do k=1,3
1457 c$$$            dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1458 c$$$          enddo
1459 c$$$          dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1460 c$$$c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1461 c$$$c    &             esclocbi,ss,ssd
1462 c$$$          escloci=ss*escloci+(1.0d0-ss)*esclocbi
1463 c$$$c         escloci=esclocbi
1464 c$$$c         write (iout,*) escloci
1465 c$$$        else if (x(2).lt.delta) then
1466 c$$$          xtemp(1)=x(1)
1467 c$$$          xtemp(2)=delta
1468 c$$$          xtemp(3)=x(3)
1469 c$$$          call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1470 c$$$          xtemp(2)=0.0d0
1471 c$$$          call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1472 c$$$          call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1473 c$$$     &        escloci,dersc(2))
1474 c$$$          call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1475 c$$$     &        ddersc0(1),dersc(1))
1476 c$$$          call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1477 c$$$     &        ddersc0(3),dersc(3))
1478 c$$$          xtemp(2)=delta
1479 c$$$          call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1480 c$$$          xtemp(2)=0.0d0
1481 c$$$          call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1482 c$$$          call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1483 c$$$     &            dersc0(2),esclocbi,dersc02)
1484 c$$$          call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1485 c$$$     &            dersc12,dersc01)
1486 c$$$          dersc0(1)=dersc01
1487 c$$$          dersc0(2)=dersc02
1488 c$$$          dersc0(3)=0.0d0
1489 c$$$          call splinthet(x(2),0.5d0*delta,ss,ssd)
1490 c$$$          do k=1,3
1491 c$$$            dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1492 c$$$          enddo
1493 c$$$          dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1494 c$$$c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1495 c$$$c    &             esclocbi,ss,ssd
1496 c$$$          escloci=ss*escloci+(1.0d0-ss)*esclocbi
1497 c$$$c         write (iout,*) escloci
1498 c$$$        else
1499 c$$$          call enesc(x,escloci,dersc,ddummy,.false.)
1500 c$$$        endif
1501 c$$$
1502 c$$$        escloc=escloc+escloci
1503 c$$$        if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1504 c$$$     &     'escloc',i,escloci
1505 c$$$c       write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1506 c$$$
1507 c$$$        gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1508 c$$$     &   wscloc*dersc(1)
1509 c$$$        gloc(ialph(i,1),icg)=wscloc*dersc(2)
1510 c$$$        gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1511 c$$$    1   continue
1512 c$$$      ENDIF
1513 c$$$      enddo
1514 c$$$      return
1515 c$$$      end
1516 c$$$
1517 c$$$C-----------------------------------------------------------------------------
1518 c$$$
1519 c$$$      subroutine egb_ij(i_sc,j_sc,evdw)
1520 c$$$C
1521 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1522 c$$$C assuming the Gay-Berne potential of interaction.
1523 c$$$C
1524 c$$$      implicit real*8 (a-h,o-z)
1525 c$$$      include 'DIMENSIONS'
1526 c$$$      include 'COMMON.GEO'
1527 c$$$      include 'COMMON.VAR'
1528 c$$$      include 'COMMON.LOCAL'
1529 c$$$      include 'COMMON.CHAIN'
1530 c$$$      include 'COMMON.DERIV'
1531 c$$$      include 'COMMON.NAMES'
1532 c$$$      include 'COMMON.INTERACT'
1533 c$$$      include 'COMMON.IOUNITS'
1534 c$$$      include 'COMMON.CALC'
1535 c$$$      include 'COMMON.CONTROL'
1536 c$$$      logical lprn
1537 c$$$      evdw=0.0D0
1538 c$$$      energy_dec=.false.
1539 c$$$c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1540 c$$$      evdw=0.0D0
1541 c$$$      lprn=.false.
1542 c$$$      ind=0
1543 c$$$c$$$      do i=iatsc_s,iatsc_e
1544 c$$$      i=i_sc
1545 c$$$        itypi=itype(i)
1546 c$$$        itypi1=itype(i+1)
1547 c$$$        xi=c(1,nres+i)
1548 c$$$        yi=c(2,nres+i)
1549 c$$$        zi=c(3,nres+i)
1550 c$$$        dxi=dc_norm(1,nres+i)
1551 c$$$        dyi=dc_norm(2,nres+i)
1552 c$$$        dzi=dc_norm(3,nres+i)
1553 c$$$c        dsci_inv=dsc_inv(itypi)
1554 c$$$        dsci_inv=vbld_inv(i+nres)
1555 c$$$c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1556 c$$$c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1557 c$$$C
1558 c$$$C Calculate SC interaction energy.
1559 c$$$C
1560 c$$$c$$$        do iint=1,nint_gr(i)
1561 c$$$c$$$          do j=istart(i,iint),iend(i,iint)
1562 c$$$        j=j_sc
1563 c$$$            ind=ind+1
1564 c$$$            itypj=itype(j)
1565 c$$$c            dscj_inv=dsc_inv(itypj)
1566 c$$$            dscj_inv=vbld_inv(j+nres)
1567 c$$$c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1568 c$$$c     &       1.0d0/vbld(j+nres)
1569 c$$$c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1570 c$$$            sig0ij=sigma(itypi,itypj)
1571 c$$$            chi1=chi(itypi,itypj)
1572 c$$$            chi2=chi(itypj,itypi)
1573 c$$$            chi12=chi1*chi2
1574 c$$$            chip1=chip(itypi)
1575 c$$$            chip2=chip(itypj)
1576 c$$$            chip12=chip1*chip2
1577 c$$$            alf1=alp(itypi)
1578 c$$$            alf2=alp(itypj)
1579 c$$$            alf12=0.5D0*(alf1+alf2)
1580 c$$$C For diagnostics only!!!
1581 c$$$c           chi1=0.0D0
1582 c$$$c           chi2=0.0D0
1583 c$$$c           chi12=0.0D0
1584 c$$$c           chip1=0.0D0
1585 c$$$c           chip2=0.0D0
1586 c$$$c           chip12=0.0D0
1587 c$$$c           alf1=0.0D0
1588 c$$$c           alf2=0.0D0
1589 c$$$c           alf12=0.0D0
1590 c$$$            xj=c(1,nres+j)-xi
1591 c$$$            yj=c(2,nres+j)-yi
1592 c$$$            zj=c(3,nres+j)-zi
1593 c$$$            dxj=dc_norm(1,nres+j)
1594 c$$$            dyj=dc_norm(2,nres+j)
1595 c$$$            dzj=dc_norm(3,nres+j)
1596 c$$$c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1597 c$$$c            write (iout,*) "j",j," dc_norm",
1598 c$$$c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1599 c$$$            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1600 c$$$            rij=dsqrt(rrij)
1601 c$$$C Calculate angle-dependent terms of energy and contributions to their
1602 c$$$C derivatives.
1603 c$$$            call sc_angular
1604 c$$$            sigsq=1.0D0/sigsq
1605 c$$$            sig=sig0ij*dsqrt(sigsq)
1606 c$$$            rij_shift=1.0D0/rij-sig+sig0ij
1607 c$$$c for diagnostics; uncomment
1608 c$$$c            rij_shift=1.2*sig0ij
1609 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1610 c$$$            if (rij_shift.le.0.0D0) then
1611 c$$$              evdw=1.0D20
1612 c$$$cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1613 c$$$cd     &        restyp(itypi),i,restyp(itypj),j,
1614 c$$$cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1615 c$$$              return
1616 c$$$            endif
1617 c$$$            sigder=-sig*sigsq
1618 c$$$c---------------------------------------------------------------
1619 c$$$            rij_shift=1.0D0/rij_shift
1620 c$$$            fac=rij_shift**expon
1621 c$$$            e1=fac*fac*aa(itypi,itypj)
1622 c$$$            e2=fac*bb(itypi,itypj)
1623 c$$$            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1624 c$$$            eps2der=evdwij*eps3rt
1625 c$$$            eps3der=evdwij*eps2rt
1626 c$$$c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1627 c$$$c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1628 c$$$            evdwij=evdwij*eps2rt*eps3rt
1629 c$$$            evdw=evdw+evdwij
1630 c$$$            if (lprn) then
1631 c$$$            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1632 c$$$            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1633 c$$$            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1634 c$$$     &        restyp(itypi),i,restyp(itypj),j,
1635 c$$$     &        epsi,sigm,chi1,chi2,chip1,chip2,
1636 c$$$     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1637 c$$$     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1638 c$$$     &        evdwij
1639 c$$$            endif
1640 c$$$
1641 c$$$            if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1642 c$$$     &                        'evdw',i,j,evdwij
1643 c$$$
1644 c$$$C Calculate gradient components.
1645 c$$$            e1=e1*eps1*eps2rt**2*eps3rt**2
1646 c$$$            fac=-expon*(e1+evdwij)*rij_shift
1647 c$$$            sigder=fac*sigder
1648 c$$$            fac=rij*fac
1649 c$$$c            fac=0.0d0
1650 c$$$C Calculate the radial part of the gradient
1651 c$$$            gg(1)=xj*fac
1652 c$$$            gg(2)=yj*fac
1653 c$$$            gg(3)=zj*fac
1654 c$$$C Calculate angular part of the gradient.
1655 c$$$            call sc_grad
1656 c$$$c$$$          enddo      ! j
1657 c$$$c$$$        enddo        ! iint
1658 c$$$c$$$      enddo          ! i
1659 c$$$      energy_dec=.false.
1660 c$$$      return
1661 c$$$      end
1662 c$$$
1663 c$$$C-----------------------------------------------------------------------------
1664 c$$$
1665 c$$$      subroutine perturb_side_chain(i,angle)
1666 c$$$      implicit none
1667 c$$$
1668 c$$$c     Includes
1669 c$$$      include 'DIMENSIONS'
1670 c$$$      include 'COMMON.CHAIN'
1671 c$$$      include 'COMMON.GEO'
1672 c$$$      include 'COMMON.VAR'
1673 c$$$      include 'COMMON.LOCAL'
1674 c$$$      include 'COMMON.IOUNITS'
1675 c$$$
1676 c$$$c     External functions
1677 c$$$      external ran_number
1678 c$$$      double precision ran_number
1679 c$$$
1680 c$$$c     Input arguments
1681 c$$$      integer i
1682 c$$$      double precision angle    ! In degrees
1683 c$$$
1684 c$$$c     Local variables
1685 c$$$      integer i_sc
1686 c$$$      double precision rad_ang,rand_v(3),length,cost,sint
1687 c$$$
1688 c$$$
1689 c$$$      i_sc=i+nres
1690 c$$$      rad_ang=angle*deg2rad
1691 c$$$
1692 c$$$      length=0.0
1693 c$$$      do while (length.lt.0.01)
1694 c$$$        rand_v(1)=ran_number(0.01D0,1.0D0)
1695 c$$$        rand_v(2)=ran_number(0.01D0,1.0D0)
1696 c$$$        rand_v(3)=ran_number(0.01D0,1.0D0)
1697 c$$$        length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1698 c$$$     +       rand_v(3)*rand_v(3)
1699 c$$$        length=sqrt(length)
1700 c$$$        rand_v(1)=rand_v(1)/length
1701 c$$$        rand_v(2)=rand_v(2)/length
1702 c$$$        rand_v(3)=rand_v(3)/length
1703 c$$$        cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1704 c$$$     +       rand_v(3)*dc_norm(3,i_sc)
1705 c$$$        length=1.0D0-cost*cost
1706 c$$$        if (length.lt.0.0D0) length=0.0D0
1707 c$$$        length=sqrt(length)
1708 c$$$        rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1709 c$$$        rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1710 c$$$        rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1711 c$$$      enddo
1712 c$$$      rand_v(1)=rand_v(1)/length
1713 c$$$      rand_v(2)=rand_v(2)/length
1714 c$$$      rand_v(3)=rand_v(3)/length
1715 c$$$
1716 c$$$      cost=dcos(rad_ang)
1717 c$$$      sint=dsin(rad_ang)
1718 c$$$      dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1719 c$$$      dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1720 c$$$      dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1721 c$$$      dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1722 c$$$      dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1723 c$$$      dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1724 c$$$      c(1,i_sc)=c(1,i)+dc(1,i_sc)
1725 c$$$      c(2,i_sc)=c(2,i)+dc(2,i_sc)
1726 c$$$      c(3,i_sc)=c(3,i)+dc(3,i_sc)
1727 c$$$
1728 c$$$      call chainbuild_cart
1729 c$$$
1730 c$$$      return
1731 c$$$      end
1732 c$$$
1733 c$$$c----------------------------------------------------------------------------
1734 c$$$
1735 c$$$      subroutine ss_relax3(i_in,j_in)
1736 c$$$      implicit none
1737 c$$$
1738 c$$$c     Includes
1739 c$$$      include 'DIMENSIONS'
1740 c$$$      include 'COMMON.VAR'
1741 c$$$      include 'COMMON.CHAIN'
1742 c$$$      include 'COMMON.IOUNITS'
1743 c$$$      include 'COMMON.INTERACT'
1744 c$$$
1745 c$$$c     External functions
1746 c$$$      external ran_number
1747 c$$$      double precision ran_number
1748 c$$$
1749 c$$$c     Input arguments
1750 c$$$      integer i_in,j_in
1751 c$$$
1752 c$$$c     Local variables
1753 c$$$      double precision energy_sc(0:n_ene),etot
1754 c$$$      double precision org_dc(3),org_dc_norm(3),org_c(3)
1755 c$$$      double precision ang_pert,rand_fact,exp_fact,beta
1756 c$$$      integer n,i_pert,i
1757 c$$$      logical notdone
1758 c$$$
1759 c$$$
1760 c$$$      beta=1.0D0
1761 c$$$
1762 c$$$      mask_r=.true.
1763 c$$$      do i=nnt,nct
1764 c$$$        mask_side(i)=0
1765 c$$$      enddo
1766 c$$$      mask_side(i_in)=1
1767 c$$$      mask_side(j_in)=1
1768 c$$$
1769 c$$$      call etotal_sc(energy_sc)
1770 c$$$      etot=energy_sc(0)
1771 c$$$c      write(iout,'(a,3d15.5)')"     SS_MC_START ",energy_sc(0),
1772 c$$$c     +     energy_sc(1),energy_sc(12)
1773 c$$$
1774 c$$$      notdone=.true.
1775 c$$$      n=0
1776 c$$$      do while (notdone)
1777 c$$$        if (mod(n,2).eq.0) then
1778 c$$$          i_pert=i_in
1779 c$$$        else
1780 c$$$          i_pert=j_in
1781 c$$$        endif
1782 c$$$        n=n+1
1783 c$$$
1784 c$$$        do i=1,3
1785 c$$$          org_dc(i)=dc(i,i_pert+nres)
1786 c$$$          org_dc_norm(i)=dc_norm(i,i_pert+nres)
1787 c$$$          org_c(i)=c(i,i_pert+nres)
1788 c$$$        enddo
1789 c$$$        ang_pert=ran_number(0.0D0,3.0D0)
1790 c$$$        call perturb_side_chain(i_pert,ang_pert)
1791 c$$$        call etotal_sc(energy_sc)
1792 c$$$        exp_fact=exp(beta*(etot-energy_sc(0)))
1793 c$$$        rand_fact=ran_number(0.0D0,1.0D0)
1794 c$$$        if (rand_fact.lt.exp_fact) then
1795 c$$$c          write(iout,'(a,3d15.5)')"     SS_MC_ACCEPT ",energy_sc(0),
1796 c$$$c     +     energy_sc(1),energy_sc(12)
1797 c$$$          etot=energy_sc(0)
1798 c$$$        else
1799 c$$$c          write(iout,'(a,3d15.5)')"     SS_MC_REJECT ",energy_sc(0),
1800 c$$$c     +     energy_sc(1),energy_sc(12)
1801 c$$$          do i=1,3
1802 c$$$            dc(i,i_pert+nres)=org_dc(i)
1803 c$$$            dc_norm(i,i_pert+nres)=org_dc_norm(i)
1804 c$$$            c(i,i_pert+nres)=org_c(i)
1805 c$$$          enddo
1806 c$$$        endif
1807 c$$$
1808 c$$$        if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1809 c$$$      enddo
1810 c$$$
1811 c$$$      mask_r=.false.
1812 c$$$
1813 c$$$      return
1814 c$$$      end
1815 c$$$
1816 c$$$c----------------------------------------------------------------------------
1817 c$$$
1818 c$$$      subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1819 c$$$      implicit none
1820 c$$$      include 'DIMENSIONS'
1821 c$$$      integer liv,lv
1822 c$$$      parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1823 c$$$*********************************************************************
1824 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1825 c$$$* the calling subprogram.                                           *
1826 c$$$* when d(i)=1.0, then v(35) is the length of the initial step,      *
1827 c$$$* calculated in the usual pythagorean way.                          *
1828 c$$$* absolute convergence occurs when the function is within v(31) of  *
1829 c$$$* zero. unless you know the minimum value in advance, abs convg     *
1830 c$$$* is probably not useful.                                           *
1831 c$$$* relative convergence is when the model predicts that the function *
1832 c$$$* will decrease by less than v(32)*abs(fun).                        *
1833 c$$$*********************************************************************
1834 c$$$      include 'COMMON.IOUNITS'
1835 c$$$      include 'COMMON.VAR'
1836 c$$$      include 'COMMON.GEO'
1837 c$$$      include 'COMMON.MINIM'
1838 c$$$      include 'COMMON.CHAIN'
1839 c$$$
1840 c$$$      double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1841 c$$$      common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1842 c$$$     +     orig_ss_dist(maxres2,maxres2)
1843 c$$$
1844 c$$$      double precision etot
1845 c$$$      integer iretcode,nfun,i_in,j_in
1846 c$$$
1847 c$$$      external dist
1848 c$$$      double precision dist
1849 c$$$      external ss_func,fdum
1850 c$$$      double precision ss_func,fdum
1851 c$$$
1852 c$$$      integer iv(liv),uiparm(2)
1853 c$$$      double precision v(lv),x(maxres6),d(maxres6),rdum
1854 c$$$      integer i,j,k
1855 c$$$
1856 c$$$
1857 c$$$      call deflt(2,iv,liv,lv,v)
1858 c$$$* 12 means fresh start, dont call deflt
1859 c$$$      iv(1)=12
1860 c$$$* max num of fun calls
1861 c$$$      if (maxfun.eq.0) maxfun=500
1862 c$$$      iv(17)=maxfun
1863 c$$$* max num of iterations
1864 c$$$      if (maxmin.eq.0) maxmin=1000
1865 c$$$      iv(18)=maxmin
1866 c$$$* controls output
1867 c$$$      iv(19)=2
1868 c$$$* selects output unit
1869 c$$$c      iv(21)=iout
1870 c$$$      iv(21)=0
1871 c$$$* 1 means to print out result
1872 c$$$      iv(22)=0
1873 c$$$* 1 means to print out summary stats
1874 c$$$      iv(23)=0
1875 c$$$* 1 means to print initial x and d
1876 c$$$      iv(24)=0
1877 c$$$* min val for v(radfac) default is 0.1
1878 c$$$      v(24)=0.1D0
1879 c$$$* max val for v(radfac) default is 4.0
1880 c$$$      v(25)=2.0D0
1881 c$$$c     v(25)=4.0D0
1882 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1883 c$$$* the sumsl default is 0.1
1884 c$$$      v(26)=0.1D0
1885 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1886 c$$$* the sumsl default is 100*machep
1887 c$$$      v(34)=v(34)/100.0D0
1888 c$$$* absolute convergence
1889 c$$$      if (tolf.eq.0.0D0) tolf=1.0D-4
1890 c$$$      v(31)=tolf
1891 c$$$      v(31)=1.0D-1
1892 c$$$* relative convergence
1893 c$$$      if (rtolf.eq.0.0D0) rtolf=1.0D-4
1894 c$$$      v(32)=rtolf
1895 c$$$      v(32)=1.0D-1
1896 c$$$* controls initial step size
1897 c$$$      v(35)=1.0D-1
1898 c$$$* large vals of d correspond to small components of step
1899 c$$$      do i=1,6*nres
1900 c$$$        d(i)=1.0D0
1901 c$$$      enddo
1902 c$$$
1903 c$$$      do i=0,2*nres
1904 c$$$        do j=1,3
1905 c$$$          orig_ss_dc(j,i)=dc(j,i)
1906 c$$$        enddo
1907 c$$$      enddo
1908 c$$$      call geom_to_var(nvar,orig_ss_var)
1909 c$$$
1910 c$$$      do i=1,nres
1911 c$$$        do j=i,nres
1912 c$$$          orig_ss_dist(j,i)=dist(j,i)
1913 c$$$          orig_ss_dist(j+nres,i)=dist(j+nres,i)
1914 c$$$          orig_ss_dist(j,i+nres)=dist(j,i+nres)
1915 c$$$          orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1916 c$$$        enddo
1917 c$$$      enddo
1918 c$$$
1919 c$$$      k=0
1920 c$$$      do i=1,nres-1
1921 c$$$        do j=1,3
1922 c$$$          k=k+1
1923 c$$$          x(k)=dc(j,i)
1924 c$$$        enddo
1925 c$$$      enddo
1926 c$$$      do i=2,nres-1
1927 c$$$        if (ialph(i,1).gt.0) then
1928 c$$$        do j=1,3
1929 c$$$          k=k+1
1930 c$$$          x(k)=dc(j,i+nres)
1931 c$$$        enddo
1932 c$$$        endif
1933 c$$$      enddo
1934 c$$$
1935 c$$$      uiparm(1)=i_in
1936 c$$$      uiparm(2)=j_in
1937 c$$$      call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1938 c$$$      etot=v(10)
1939 c$$$      iretcode=iv(1)
1940 c$$$      nfun=iv(6)+iv(30)
1941 c$$$
1942 c$$$      k=0
1943 c$$$      do i=1,nres-1
1944 c$$$        do j=1,3
1945 c$$$          k=k+1
1946 c$$$          dc(j,i)=x(k)
1947 c$$$        enddo
1948 c$$$      enddo
1949 c$$$      do i=2,nres-1
1950 c$$$        if (ialph(i,1).gt.0) then
1951 c$$$        do j=1,3
1952 c$$$          k=k+1
1953 c$$$          dc(j,i+nres)=x(k)
1954 c$$$        enddo
1955 c$$$        endif
1956 c$$$      enddo
1957 c$$$      call chainbuild_cart
1958 c$$$
1959 c$$$      return
1960 c$$$      end
1961 c$$$
1962 c$$$C-----------------------------------------------------------------------------
1963 c$$$
1964 c$$$      subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1965 c$$$      implicit none
1966 c$$$      include 'DIMENSIONS'
1967 c$$$      include 'COMMON.DERIV'
1968 c$$$      include 'COMMON.IOUNITS'
1969 c$$$      include 'COMMON.VAR'
1970 c$$$      include 'COMMON.CHAIN'
1971 c$$$      include 'COMMON.INTERACT'
1972 c$$$      include 'COMMON.SBRIDGE'
1973 c$$$
1974 c$$$      double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1975 c$$$      common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1976 c$$$     +     orig_ss_dist(maxres2,maxres2)
1977 c$$$
1978 c$$$      integer n
1979 c$$$      double precision x(maxres6)
1980 c$$$      integer nf
1981 c$$$      double precision f
1982 c$$$      integer uiparm(2)
1983 c$$$      real*8 urparm(1)
1984 c$$$      external ufparm
1985 c$$$      double precision ufparm
1986 c$$$
1987 c$$$      external dist
1988 c$$$      double precision dist
1989 c$$$
1990 c$$$      integer i,j,k,ss_i,ss_j
1991 c$$$      double precision tempf,var(maxvar)
1992 c$$$
1993 c$$$
1994 c$$$      ss_i=uiparm(1)
1995 c$$$      ss_j=uiparm(2)
1996 c$$$      f=0.0D0
1997 c$$$
1998 c$$$      k=0
1999 c$$$      do i=1,nres-1
2000 c$$$        do j=1,3
2001 c$$$          k=k+1
2002 c$$$          dc(j,i)=x(k)
2003 c$$$        enddo
2004 c$$$      enddo
2005 c$$$      do i=2,nres-1
2006 c$$$        if (ialph(i,1).gt.0) then
2007 c$$$        do j=1,3
2008 c$$$          k=k+1
2009 c$$$          dc(j,i+nres)=x(k)
2010 c$$$        enddo
2011 c$$$        endif
2012 c$$$      enddo
2013 c$$$      call chainbuild_cart
2014 c$$$
2015 c$$$      call geom_to_var(nvar,var)
2016 c$$$
2017 c$$$c     Constraints on all angles
2018 c$$$      do i=1,nvar
2019 c$$$        tempf=var(i)-orig_ss_var(i)
2020 c$$$        f=f+tempf*tempf
2021 c$$$      enddo
2022 c$$$
2023 c$$$c     Constraints on all distances
2024 c$$$      do i=1,nres-1
2025 c$$$        if (i.gt.1) then
2026 c$$$          tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
2027 c$$$          f=f+tempf*tempf
2028 c$$$        endif
2029 c$$$        do j=i+1,nres
2030 c$$$          tempf=dist(j,i)-orig_ss_dist(j,i)
2031 c$$$          if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
2032 c$$$          tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
2033 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
2034 c$$$          tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
2035 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
2036 c$$$          tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
2037 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
2038 c$$$        enddo
2039 c$$$      enddo
2040 c$$$
2041 c$$$c     Constraints for the relevant CYS-CYS
2042 c$$$      tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
2043 c$$$      f=f+tempf*tempf
2044 c$$$CCCCCCCCCCCCCCCCC      ADD SOME ANGULAR STUFF
2045 c$$$
2046 c$$$c$$$      if (nf.ne.nfl) then
2047 c$$$c$$$        write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
2048 c$$$c$$$     +       f,dist(5+nres,14+nres)
2049 c$$$c$$$      endif
2050 c$$$
2051 c$$$      nfl=nf
2052 c$$$
2053 c$$$      return
2054 c$$$      end
2055 c$$$
2056 c$$$C-----------------------------------------------------------------------------
2057 c$$$C-----------------------------------------------------------------------------
2058          subroutine triple_ssbond_ene(resi,resj,resk,eij)
2059       include 'DIMENSIONS'
2060       include 'COMMON.SBRIDGE'
2061       include 'COMMON.CHAIN'
2062       include 'COMMON.DERIV'
2063       include 'COMMON.LOCAL'
2064       include 'COMMON.INTERACT'
2065       include 'COMMON.VAR'
2066       include 'COMMON.IOUNITS'
2067       include 'COMMON.CALC'
2068 #ifndef CLUST
2069 #ifndef WHAM
2070 C      include 'COMMON.MD'
2071 #endif
2072 #endif
2073
2074 c     External functions
2075       double precision h_base
2076       external h_base
2077
2078 c     Input arguments
2079       integer resi,resj,resk
2080
2081 c     Output arguments
2082       double precision eij,eij1,eij2,eij3
2083
2084 c     Local variables
2085       logical havebond
2086 c      integer itypi,itypj,k,l
2087       double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2088       double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2089       double precision xik,yik,zik,xjk,yjk,zjk
2090       double precision sig0ij,ljd,sig,fac,e1,e2
2091       double precision dcosom1(3),dcosom2(3),ed
2092       double precision pom1,pom2
2093       double precision ljA,ljB,ljXs
2094       double precision d_ljB(1:3)
2095       double precision ssA,ssB,ssC,ssXs
2096       double precision ssxm,ljxm,ssm,ljm
2097       double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2098
2099       i=resi
2100       j=resj
2101       k=resk
2102 C      write(iout,*) resi,resj,resk
2103       itypi=itype(i)
2104       dxi=dc_norm(1,nres+i)
2105       dyi=dc_norm(2,nres+i)
2106       dzi=dc_norm(3,nres+i)
2107       dsci_inv=vbld_inv(i+nres)
2108       xi=c(1,nres+i)
2109       yi=c(2,nres+i)
2110       zi=c(3,nres+i)
2111
2112       itypj=itype(j)
2113       xj=c(1,nres+j)
2114       yj=c(2,nres+j)
2115       zj=c(3,nres+j)
2116
2117       dxj=dc_norm(1,nres+j)
2118       dyj=dc_norm(2,nres+j)
2119       dzj=dc_norm(3,nres+j)
2120       dscj_inv=vbld_inv(j+nres)
2121       itypk=itype(k)
2122       xk=c(1,nres+k)
2123       yk=c(2,nres+k)
2124       zk=c(3,nres+k)
2125
2126       dxk=dc_norm(1,nres+k)
2127       dyk=dc_norm(2,nres+k)
2128       dzk=dc_norm(3,nres+k)
2129       dscj_inv=vbld_inv(k+nres)
2130       xij=xj-xi
2131       xik=xk-xi
2132       xjk=xk-xj
2133       yij=yj-yi
2134       yik=yk-yi
2135       yjk=yk-yj
2136       zij=zj-zi
2137       zik=zk-zi
2138       zjk=zk-zj
2139       rrij=(xij*xij+yij*yij+zij*zij)
2140       rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
2141       rrik=(xik*xik+yik*yik+zik*zik)
2142       rik=dsqrt(rrik)
2143       rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2144       rjk=dsqrt(rrjk)
2145 C there are three combination of distances for each trisulfide bonds
2146 C The first case the ith atom is the center
2147 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2148 C distance y is second distance the a,b,c,d are parameters derived for
2149 C this problem d parameter was set as a penalty currenlty set to 1.
2150       eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2151 C second case jth atom is center
2152       eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2153 C the third case kth atom is the center
2154       eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2155 C      eij2=0.0
2156 C      eij3=0.0
2157 C      eij1=0.0
2158       eij=eij1+eij2+eij3
2159 C      write(iout,*)i,j,k,eij
2160 C The energy penalty calculated now time for the gradient part
2161 C derivative over rij
2162       fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2163      &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2164             gg(1)=xij*fac/rij
2165             gg(2)=yij*fac/rij
2166             gg(3)=zij*fac/rij
2167       do m=1,3
2168         gvdwx(m,i)=gvdwx(m,i)-gg(m)
2169         gvdwx(m,j)=gvdwx(m,j)+gg(m)
2170       enddo
2171       do l=1,3
2172         gvdwc(l,i)=gvdwc(l,i)-gg(l)
2173         gvdwc(l,j)=gvdwc(l,j)+gg(l)
2174       enddo
2175 C now derivative over rik
2176       fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2177      &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2178             gg(1)=xik*fac/rik
2179             gg(2)=yik*fac/rik
2180             gg(3)=zik*fac/rik
2181       do m=1,3
2182         gvdwx(m,i)=gvdwx(m,i)-gg(m)
2183         gvdwx(m,k)=gvdwx(m,k)+gg(m)
2184       enddo
2185       do l=1,3
2186         gvdwc(l,i)=gvdwc(l,i)-gg(l)
2187         gvdwc(l,k)=gvdwc(l,k)+gg(l)
2188       enddo
2189 C now derivative over rjk
2190       fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2191      &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2192             gg(1)=xjk*fac/rjk
2193             gg(2)=yjk*fac/rjk
2194             gg(3)=zjk*fac/rjk
2195       do m=1,3
2196         gvdwx(m,j)=gvdwx(m,j)-gg(m)
2197         gvdwx(m,k)=gvdwx(m,k)+gg(m)
2198       enddo
2199       do l=1,3
2200         gvdwc(l,j)=gvdwc(l,j)-gg(l)
2201         gvdwc(l,k)=gvdwc(l,k)+gg(l)
2202       enddo
2203       return
2204       end