source/maxlik/src_MD_T_maxlik-NEWCORR.safe/ssMD.F

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