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