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