source/wham/src-M/ssMD.F

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