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