source/unres/src_MD-M/ssMD.F

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