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