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