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