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

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