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