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