source/wham/src-HCD-5D/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       subroutine dyn_ssbond_ene(resi,resj,eij)
  86       implicit none
  87       include 'DIMENSIONS'
  88       include 'DIMENSIONS.ZSCOPT'
  89       include 'COMMON.SBRIDGE'
  90       include 'COMMON.CHAIN'
  91       include 'COMMON.DERIV'
  92       include 'COMMON.LOCAL'
  93       include 'COMMON.INTERACT'
  94       include 'COMMON.VAR'
  95       include 'COMMON.IOUNITS'
  96       include 'COMMON.CALC'
  97       include 'COMMON.NAMES'
  98 #ifndef CLUST
  99 #ifndef WHAM
 100       include 'COMMON.MD'
 101 #endif
 102 #endif
 103
 104 c     External functions
 105       double precision h_base
 106       external h_base
 107
 108 c     Input arguments
 109       integer resi,resj
 110
 111 c     Output arguments
 112       double precision eij
 113
 114 c     Local variables
 115       logical havebond
 116 c      integer itypi,itypj,k,l
 117       double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
 118       double precision sig0ij,ljd,sig,fac,e1,e2
 119       double precision dcosom1(3),dcosom2(3),ed
 120       double precision pom1,pom2
 121       double precision ljA,ljB,ljXs
 122       double precision d_ljB(1:3)
 123       double precision ssA,ssB,ssC,ssXs
 124       double precision ssxm,ljxm,ssm,ljm
 125       double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
 126       double precision f1,f2,h1,h2,hd1,hd2
 127       double precision omega,delta_inv,deltasq_inv,fac1,fac2
 128 c-------FIRST METHOD
 129       double precision xm,d_xm(1:3)
 130       double precision sslipi,sslipj,ssgradlipi,ssgradlipj
 131       integer ici,icj,itypi,itypj
 132       double precision boxshift,sscale,sscagrad
 133       double precision aa,bb
 134 c-------END FIRST METHOD
 135 c-------SECOND METHOD
 136 c$$$      double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
 137 c-------END SECOND METHOD
 138
 139 c-------TESTING CODE
 140       logical checkstop,transgrad
 141       common /sschecks/ checkstop,transgrad
 142
 143       integer icheck,nicheck,jcheck,njcheck
 144       double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
 145 c-------END TESTING CODE
 146
 147
 148       i=resi
 149       j=resj
 150       ici=icys(i)
 151       icj=icys(j)
 152       if (ici.eq.0 .or. icj.eq.0) then
 153         write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
 154      &  "Attempt to create",
 155      &  " a disulfide link between non-cysteine residues ",restyp(i),i,
 156      &  restyp(j),j
 157         stop
 158       endif
 159       itypi=itype(i)
 160       dxi=dc_norm(1,nres+i)
 161       dyi=dc_norm(2,nres+i)
 162       dzi=dc_norm(3,nres+i)
 163       dsci_inv=vbld_inv(i+nres)
 164       xi=c(1,nres+i)
 165       yi=c(2,nres+i)
 166       zi=c(3,nres+i)
 167       call to_box(xi,yi,zi)
 168 C define scaling factor for lipids
 169
 170 C        if (positi.le.0) positi=positi+boxzsize
 171 C        print *,i
 172 C first for peptide groups
 173 c for each residue check if it is in lipid or lipid water border area
 174       call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
 175       itypj=itype(j)
 176       xj=c(1,nres+j)
 177       yj=c(2,nres+j)
 178       zj=c(3,nres+j)
 179       call to_box(xj,yj,zj)
 180       call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
 181       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 182      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
 183       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 184      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
 185       xj=boxshift(xj-xi,boxxsize)
 186       yj=boxshift(yj-yi,boxysize)
 187       zj=boxshift(zj-zi,boxzsize)
 188       dxj=dc_norm(1,nres+j)
 189       dyj=dc_norm(2,nres+j)
 190       dzj=dc_norm(3,nres+j)
 191       dscj_inv=vbld_inv(j+nres)
 192
 193       chi1=chi(itypi,itypj)
 194       chi2=chi(itypj,itypi)
 195       chi12=chi1*chi2
 196       chip1=chip(itypi)
 197       chip2=chip(itypj)
 198       chip12=chip1*chip2
 199       alf1=alp(itypi)
 200       alf2=alp(itypj)
 201       alf12=0.5D0*(alf1+alf2)
 202
 203       rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 204       rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
 205             sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
 206             sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
 207 c     The following are set in sc_angular
 208 c      erij(1)=xj*rij
 209 c      erij(2)=yj*rij
 210 c      erij(3)=zj*rij
 211 c      om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
 212 c      om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
 213 c      om12=dxi*dxj+dyi*dyj+dzi*dzj
 214       call sc_angular
 215       rij=1.0D0/rij  ! Reset this so it makes sense
 216
 217       sig0ij=sigma(itypi,itypj)
 218       sig=sig0ij*dsqrt(1.0D0/sigsq)
 219
 220       ljXs=sig-sig0ij
 221       ljA=eps1*eps2rt**2*eps3rt**2
 222       ljB=ljA*bb
 223       ljA=ljA*aa
 224       ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
 225
 226       ssXs=d0cm
 227       deltat1=1.0d0-om1
 228       deltat2=1.0d0+om2
 229       deltat12=om2-om1+2.0d0
 230       cosphi=om12-om1*om2
 231       ssA=akcm
 232       ssB=akct*deltat12
 233       ssC=ss_depth
 234      &     +akth*(deltat1*deltat1+deltat2*deltat2)
 235      &     +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
 236       ssxm=ssXs-0.5D0*ssB/ssA
 237
 238 c-------TESTING CODE
 239 c$$$c     Some extra output
 240 c$$$      ssm=ssC-0.25D0*ssB*ssB/ssA
 241 c$$$      ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
 242 c$$$      ssx0=ssB*ssB-4.0d0*ssA*ssC
 243 c$$$      if (ssx0.gt.0.0d0) then
 244 c$$$        ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
 245 c$$$      else
 246 c$$$        ssx0=ssxm
 247 c$$$      endif
 248 c$$$      ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 249 c$$$      write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
 250 c$$$     &     ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
 251 c$$$      return
 252 c-------END TESTING CODE
 253
 254 c-------TESTING CODE
 255 c     Stop and plot energy and derivative as a function of distance
 256       if (checkstop) then
 257         ssm=ssC-0.25D0*ssB*ssB/ssA
 258         ljm=-0.25D0*ljB*bb/aa
 259         if (ssm.lt.ljm .and.
 260      &       dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
 261           nicheck=1000
 262           njcheck=1
 263           deps=0.5d-7
 264         else
 265           checkstop=.false.
 266         endif
 267       endif
 268       if (.not.checkstop) then
 269         nicheck=0
 270         njcheck=-1
 271       endif
 272
 273       do icheck=0,nicheck
 274       do jcheck=-1,njcheck
 275       if (checkstop) rij=(ssxm-1.0d0)+
 276      &       ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
 277 c-------END TESTING CODE
 278
 279       if (rij.gt.ljxm) then
 280         havebond=.false.
 281         ljd=rij-ljXs
 282         fac=(1.0D0/ljd)**expon
 283         e1=fac*fac*aa
 284         e2=fac*bb
 285         eij=eps1*eps2rt*eps3rt*(e1+e2)
 286         eps2der=eij*eps3rt
 287         eps3der=eij*eps2rt
 288         eij=eij*eps2rt*eps3rt*sss
 289
 290         sigder=-sig/sigsq
 291         e1=e1*eps1*eps2rt**2*eps3rt**2
 292         ed=-expon*(e1+eij)/ljd
 293         sigder=ed*sigder
 294         ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
 295         eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
 296         eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
 297         eom12=eij*eps1_om12+eps2der*eps2rt_om12
 298      &       -2.0D0*alf12*eps3der+sigder*sigsq_om12
 299       else if (rij.lt.ssxm) then
 300         havebond=.true.
 301         ssd=rij-ssXs
 302         eij=ssA*ssd*ssd+ssB*ssd+ssC
 303         eij=eij*sss
 304         ed=2*akcm*ssd+akct*deltat12
 305         ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
 306         pom1=akct*ssd
 307         pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
 308         eom1=-2*akth*deltat1-pom1-om2*pom2
 309         eom2= 2*akth*deltat2+pom1-om1*pom2
 310         eom12=pom2
 311       else
 312         omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
 313
 314         d_ssxm(1)=0.5D0*akct/ssA
 315         d_ssxm(2)=-d_ssxm(1)
 316         d_ssxm(3)=0.0D0
 317
 318         d_ljxm(1)=sig0ij/sqrt(sigsq**3)
 319         d_ljxm(2)=d_ljxm(1)*sigsq_om2
 320         d_ljxm(3)=d_ljxm(1)*sigsq_om12
 321         d_ljxm(1)=d_ljxm(1)*sigsq_om1
 322
 323 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
 324         xm=0.5d0*(ssxm+ljxm)
 325         do k=1,3
 326           d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
 327         enddo
 328         if (rij.lt.xm) then
 329           havebond=.true.
 330           ssm=ssC-0.25D0*ssB*ssB/ssA
 331           d_ssm(1)=0.5D0*akct*ssB/ssA
 332           d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
 333           d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
 334           d_ssm(3)=omega
 335           f1=(rij-xm)/(ssxm-xm)
 336           f2=(rij-ssxm)/(xm-ssxm)
 337           h1=h_base(f1,hd1)
 338           h2=h_base(f2,hd2)
 339           eij=ssm*h1+Ht*h2
 340           delta_inv=1.0d0/(xm-ssxm)
 341           deltasq_inv=delta_inv*delta_inv
 342           fac=ssm*hd1-Ht*hd2
 343           fac1=deltasq_inv*fac*(xm-rij)
 344           fac2=deltasq_inv*fac*(rij-ssxm)
 345           ed=delta_inv*(Ht*hd2-ssm*hd1)
 346           eij=eij*sss
 347           ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
 348           eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
 349           eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
 350           eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
 351         else
 352           havebond=.false.
 353           ljm=-0.25D0*ljB*bb/aa
 354           d_ljm(1)=-0.5D0*bb/aa*ljB
 355           d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
 356           d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
 357      +         alf12/eps3rt)
 358           d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
 359           f1=(rij-ljxm)/(xm-ljxm)
 360           f2=(rij-xm)/(ljxm-xm)
 361           h1=h_base(f1,hd1)
 362           h2=h_base(f2,hd2)
 363           eij=Ht*h1+ljm*h2
 364           delta_inv=1.0d0/(ljxm-xm)
 365           deltasq_inv=delta_inv*delta_inv
 366           fac=Ht*hd1-ljm*hd2
 367           fac1=deltasq_inv*fac*(ljxm-rij)
 368           fac2=deltasq_inv*fac*(rij-xm)
 369           ed=delta_inv*(ljm*hd2-Ht*hd1)
 370           eij=eij*sss
 371           ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
 372           eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
 373           eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
 374           eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
 375         endif
 376 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
 377
 378 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 379 c$$$        ssd=rij-ssXs
 380 c$$$        ljd=rij-ljXs
 381 c$$$        fac1=rij-ljxm
 382 c$$$        fac2=rij-ssxm
 383 c$$$
 384 c$$$        d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
 385 c$$$        d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
 386 c$$$        d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
 387 c$$$
 388 c$$$        ssm=ssC-0.25D0*ssB*ssB/ssA
 389 c$$$        d_ssm(1)=0.5D0*akct*ssB/ssA
 390 c$$$        d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
 391 c$$$        d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
 392 c$$$        d_ssm(3)=omega
 393 c$$$
 394 c$$$        ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
 395 c$$$        do k=1,3
 396 c$$$          d_ljm(k)=ljm*d_ljB(k)
 397 c$$$        enddo
 398 c$$$        ljm=ljm*ljB
 399 c$$$
 400 c$$$        ss=ssA*ssd*ssd+ssB*ssd+ssC
 401 c$$$        d_ss(0)=2.0d0*ssA*ssd+ssB
 402 c$$$        d_ss(2)=akct*ssd
 403 c$$$        d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
 404 c$$$        d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
 405 c$$$        d_ss(3)=omega
 406 c$$$
 407 c$$$        ljf=bb(itypi,itypj)/aa(itypi,itypj)
 408 c$$$        ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
 409 c$$$        d_ljf(0)=ljf*2.0d0*ljB*fac1
 410 c$$$        do k=1,3
 411 c$$$          d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
 412 c$$$     &         2.0d0*ljB*fac1*d_ljxm(k))
 413 c$$$        enddo
 414 c$$$        ljf=ljm+ljf*ljB*fac1*fac1
 415 c$$$
 416 c$$$        f1=(rij-ljxm)/(ssxm-ljxm)
 417 c$$$        f2=(rij-ssxm)/(ljxm-ssxm)
 418 c$$$        h1=h_base(f1,hd1)
 419 c$$$        h2=h_base(f2,hd2)
 420 c$$$        eij=ss*h1+ljf*h2
 421 c$$$        delta_inv=1.0d0/(ljxm-ssxm)
 422 c$$$        deltasq_inv=delta_inv*delta_inv
 423 c$$$        fac=ljf*hd2-ss*hd1
 424 c$$$        ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
 425 c$$$        eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
 426 c$$$     &       (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
 427 c$$$        eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
 428 c$$$     &       (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
 429 c$$$        eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
 430 c$$$     &       (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
 431 c$$$
 432 c$$$        havebond=.false.
 433 c$$$        if (ed.gt.0.0d0) havebond=.true.
 434 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 435
 436       endif
 437
 438       if (havebond) then
 439 #ifndef CLUST
 440 #ifndef WHAM
 441 c        if (dyn_ssbond_ij(i,j).eq.1.0d300) then
 442 c          write(iout,'(a15,f12.2,f8.1,2i5)')
 443 c     &         "SSBOND_E_FORM",totT,t_bath,i,j
 444 c        endif
 445 #endif
 446 #endif
 447         dyn_ssbond_ij(ici,icj)=eij
 448       else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300)
 449      &then
 450         dyn_ssbond_ij(ici,icj)=1.0d300
 451 #ifndef CLUST
 452 #ifndef WHAM
 453 c        write(iout,'(a15,f12.2,f8.1,2i5)')
 454 c     &       "SSBOND_E_BREAK",totT,t_bath,i,j
 455 #endif
 456 #endif
 457       endif
 458
 459 c-------TESTING CODE
 460       if (checkstop) then
 461         if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
 462      &       "CHECKSTOP",rij,eij,ed
 463         echeck(jcheck)=eij
 464       endif
 465       enddo
 466       if (checkstop) then
 467         write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
 468       endif
 469       enddo
 470       if (checkstop) then
 471         transgrad=.true.
 472         checkstop=.false.
 473       endif
 474 c-------END TESTING CODE
 475             gg_lipi(3)=ssgradlipi*eij
 476             gg_lipj(3)=ssgradlipj*eij
 477
 478       do k=1,3
 479         dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
 480         dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
 481       enddo
 482       do k=1,3
 483         gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
 484       enddo
 485       do k=1,3
 486         gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
 487      &       +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
 488      &       +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 489         gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
 490      &       +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
 491      &       +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
 492       enddo
 493 cgrad      do k=i,j-1
 494 cgrad        do l=1,3
 495 cgrad          gvdwc(l,k)=gvdwc(l,k)+gg(l)
 496 cgrad        enddo
 497 cgrad      enddo
 498
 499       do l=1,3
 500         gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
 501         gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
 502       enddo
 503
 504       return
 505       end
 506 C-----------------------------------------------------------------------------
 507
 508       double precision function h_base(x,deriv)
 509 c     A smooth function going 0->1 in range [0,1]
 510 c     It should NOT be called outside range [0,1], it will not work there.
 511       implicit none
 512
 513 c     Input arguments
 514       double precision x
 515
 516 c     Output arguments
 517       double precision deriv
 518
 519 c     Local variables
 520       double precision xsq
 521
 522
 523 c     Two parabolas put together.  First derivative zero at extrema
 524 c$$$      if (x.lt.0.5D0) then
 525 c$$$        h_base=2.0D0*x*x
 526 c$$$        deriv=4.0D0*x
 527 c$$$      else
 528 c$$$        deriv=1.0D0-x
 529 c$$$        h_base=1.0D0-2.0D0*deriv*deriv
 530 c$$$        deriv=4.0D0*deriv
 531 c$$$      endif
 532
 533 c     Third degree polynomial.  First derivative zero at extrema
 534       h_base=x*x*(3.0d0-2.0d0*x)
 535       deriv=6.0d0*x*(1.0d0-x)
 536
 537 c     Fifth degree polynomial.  First and second derivatives zero at extrema
 538 c$$$      xsq=x*x
 539 c$$$      h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
 540 c$$$      deriv=x-1.0d0
 541 c$$$      deriv=deriv*deriv
 542 c$$$      deriv=30.0d0*xsq*deriv
 543
 544       return
 545       end
 546 c----------------------------------------------------------------------------
 547       subroutine dyn_set_nss
 548 c     Adjust nss and other relevant variables based on dyn_ssbond_ij
 549 c      implicit none
 550
 551 c     Includes
 552       include 'DIMENSIONS'
 553 #ifdef MPI
 554       include "mpif.h"
 555 #endif
 556       include 'COMMON.SBRIDGE'
 557       include 'COMMON.CHAIN'
 558       include 'COMMON.IOUNITS'
 559 C      include 'COMMON.SETUP'
 560 #ifndef CLUST
 561 #ifndef WHAM
 562 C      include 'COMMON.MD'
 563 #endif
 564 #endif
 565
 566 c     Local variables
 567       double precision emin
 568       integer i,j,imin
 569       integer diff,allflag(maxdim_cont),allnss,
 570      &     allihpb(maxdim_cont),alljhpb(maxdim_cont),
 571      &     newnss,newihpb(maxdim_cont),newjhpb(maxdim_cont)
 572       logical found
 573       integer i_newnss(1024),displ(0:1024)
 574       integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss
 575
 576       allnss=0
 577       do i=1,ns-1
 578         do j=i+1,ns
 579           if (dyn_ssbond_ij(i,j).lt.1.0d300) then
 580             allnss=allnss+1
 581             allflag(allnss)=0
 582             allihpb(allnss)=i
 583             alljhpb(allnss)=j
 584           endif
 585         enddo
 586       enddo
 587
 588 cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
 589
 590  1    emin=1.0d300
 591       do i=1,allnss
 592         if (allflag(i).eq.0 .and.
 593      &       dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
 594           emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
 595           imin=i
 596         endif
 597       enddo
 598       if (emin.lt.1.0d300) then
 599         allflag(imin)=1
 600         do i=1,allnss
 601           if (allflag(i).eq.0 .and.
 602      &         (allihpb(i).eq.allihpb(imin) .or.
 603      &         alljhpb(i).eq.allihpb(imin) .or.
 604      &         allihpb(i).eq.alljhpb(imin) .or.
 605      &         alljhpb(i).eq.alljhpb(imin))) then
 606             allflag(i)=-1
 607           endif
 608         enddo
 609         goto 1
 610       endif
 611
 612 cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
 613
 614       newnss=0
 615       do i=1,allnss
 616         if (allflag(i).eq.1) then
 617           newnss=newnss+1
 618           newihpb(newnss)=allihpb(i)
 619           newjhpb(newnss)=alljhpb(i)
 620         endif
 621       enddo
 622
 623 #ifdef MPI
 624       if (nfgtasks.gt.1)then
 625
 626         call MPI_Reduce(newnss,g_newnss,1,
 627      &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
 628         call MPI_Gather(newnss,1,MPI_INTEGER,
 629      &                  i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
 630         displ(0)=0
 631         do i=1,nfgtasks-1,1
 632           displ(i)=i_newnss(i-1)+displ(i-1)
 633         enddo
 634         call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
 635      &                   g_newihpb,i_newnss,displ,MPI_INTEGER,
 636      &                   king,FG_COMM,IERR)
 637         call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
 638      &                   g_newjhpb,i_newnss,displ,MPI_INTEGER,
 639      &                   king,FG_COMM,IERR)
 640         if(fg_rank.eq.0) then
 641 c         print *,'g_newnss',g_newnss
 642 c         print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
 643 c         print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
 644          newnss=g_newnss
 645          do i=1,newnss
 646           newihpb(i)=g_newihpb(i)
 647           newjhpb(i)=g_newjhpb(i)
 648          enddo
 649         endif
 650       endif
 651 #endif
 652
 653       diff=newnss-nss
 654
 655 cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
 656
 657       do i=1,nss
 658         found=.false.
 659         do j=1,newnss
 660           if (idssb(i).eq.newihpb(j) .and.
 661      &         jdssb(i).eq.newjhpb(j)) found=.true.
 662         enddo
 663 #ifndef CLUST
 664 #ifndef WHAM
 665 c        if (.not.found.and.fg_rank.eq.0)
 666 c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
 667 c     &       "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
 668 #endif
 669 #endif
 670       enddo
 671
 672       do i=1,newnss
 673         found=.false.
 674         do j=1,nss
 675           if (newihpb(i).eq.idssb(j) .and.
 676      &         newjhpb(i).eq.jdssb(j)) found=.true.
 677         enddo
 678 #ifndef CLUST
 679 #ifndef WHAM
 680 c        if (.not.found.and.fg_rank.eq.0)
 681 c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
 682 c     &       "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
 683 #endif
 684 #endif
 685       enddo
 686
 687       nss=newnss
 688       do i=1,nss
 689         idssb(i)=newihpb(i)
 690         jdssb(i)=newjhpb(i)
 691       enddo
 692
 693       return
 694       end
 695
 696 c----------------------------------------------------------------------------
 697
 698 #ifdef SSREAD
 699 #ifdef WHAM
 700       subroutine read_ssHist
 701       implicit none
 702
 703 c     Includes
 704       include 'DIMENSIONS'
 705       include "DIMENSIONS.FREE"
 706       include 'COMMON.FREE'
 707
 708 c     Local variables
 709       integer i,j
 710       character*80 controlcard
 711
 712       do i=1,dyn_nssHist
 713         call card_concat(controlcard,.true.)
 714         read(controlcard,*)
 715      &       dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
 716       enddo
 717
 718       return
 719       end
 720 #endif
 721 #endif
 722 c$$$C----------------------------------------------------------------------------
 723       subroutine triple_ssbond_ene(resi,resj,resk,eij)
 724       include 'DIMENSIONS'
 725       include 'COMMON.SBRIDGE'
 726       include 'COMMON.CHAIN'
 727       include 'COMMON.DERIV'
 728       include 'COMMON.LOCAL'
 729       include 'COMMON.INTERACT'
 730       include 'COMMON.VAR'
 731       include 'COMMON.IOUNITS'
 732       include 'COMMON.CALC'
 733 #ifndef CLUST
 734 #ifndef WHAM
 735 C      include 'COMMON.MD'
 736 #endif
 737 #endif
 738
 739 c     External functions
 740       double precision h_base
 741       external h_base
 742
 743 c     Input arguments
 744       integer resi,resj,resk
 745
 746 c     Output arguments
 747       double precision eij,eij1,eij2,eij3
 748
 749 c     Local variables
 750       logical havebond
 751 c      integer itypi,itypj,k,l
 752       double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
 753       double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
 754       double precision xik,yik,zik,xjk,yjk,zjk
 755       double precision sig0ij,ljd,sig,fac,e1,e2
 756       double precision dcosom1(3),dcosom2(3),ed
 757       double precision pom1,pom2
 758       double precision ljA,ljB,ljXs
 759       double precision d_ljB(1:3)
 760       double precision ssA,ssB,ssC,ssXs
 761       double precision ssxm,ljxm,ssm,ljm
 762       double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
 763
 764       i=resi
 765       j=resj
 766       k=resk
 767 C      write(iout,*) resi,resj,resk
 768       itypi=itype(i)
 769       dxi=dc_norm(1,nres+i)
 770       dyi=dc_norm(2,nres+i)
 771       dzi=dc_norm(3,nres+i)
 772       dsci_inv=vbld_inv(i+nres)
 773       xi=c(1,nres+i)
 774       yi=c(2,nres+i)
 775       zi=c(3,nres+i)
 776
 777       itypj=itype(j)
 778       xj=c(1,nres+j)
 779       yj=c(2,nres+j)
 780       zj=c(3,nres+j)
 781
 782       dxj=dc_norm(1,nres+j)
 783       dyj=dc_norm(2,nres+j)
 784       dzj=dc_norm(3,nres+j)
 785       dscj_inv=vbld_inv(j+nres)
 786       itypk=itype(k)
 787       xk=c(1,nres+k)
 788       yk=c(2,nres+k)
 789       zk=c(3,nres+k)
 790
 791       dxk=dc_norm(1,nres+k)
 792       dyk=dc_norm(2,nres+k)
 793       dzk=dc_norm(3,nres+k)
 794       dscj_inv=vbld_inv(k+nres)
 795       xij=xj-xi
 796       xik=xk-xi
 797       xjk=xk-xj
 798       yij=yj-yi
 799       yik=yk-yi
 800       yjk=yk-yj
 801       zij=zj-zi
 802       zik=zk-zi
 803       zjk=zk-zj
 804       rrij=(xij*xij+yij*yij+zij*zij)
 805       rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
 806       rrik=(xik*xik+yik*yik+zik*zik)
 807       rik=dsqrt(rrik)
 808       rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
 809       rjk=dsqrt(rrjk)
 810 C there are three combination of distances for each trisulfide bonds
 811 C The first case the ith atom is the center
 812 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
 813 C distance y is second distance the a,b,c,d are parameters derived for
 814 C this problem d parameter was set as a penalty currenlty set to 1.
 815       eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
 816 C second case jth atom is center
 817       eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
 818 C the third case kth atom is the center
 819       eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
 820 C      eij2=0.0
 821 C      eij3=0.0
 822 C      eij1=0.0
 823       eij=eij1+eij2+eij3
 824 C      write(iout,*)i,j,k,eij
 825 C The energy penalty calculated now time for the gradient part
 826 C derivative over rij
 827       fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
 828      &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
 829             gg(1)=xij*fac/rij
 830             gg(2)=yij*fac/rij
 831             gg(3)=zij*fac/rij
 832       do m=1,3
 833         gvdwx(m,i)=gvdwx(m,i)-gg(m)
 834         gvdwx(m,j)=gvdwx(m,j)+gg(m)
 835       enddo
 836       do l=1,3
 837         gvdwc(l,i)=gvdwc(l,i)-gg(l)
 838         gvdwc(l,j)=gvdwc(l,j)+gg(l)
 839       enddo
 840 C now derivative over rik
 841       fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
 842      &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
 843             gg(1)=xik*fac/rik
 844             gg(2)=yik*fac/rik
 845             gg(3)=zik*fac/rik
 846       do m=1,3
 847         gvdwx(m,i)=gvdwx(m,i)-gg(m)
 848         gvdwx(m,k)=gvdwx(m,k)+gg(m)
 849       enddo
 850       do l=1,3
 851         gvdwc(l,i)=gvdwc(l,i)-gg(l)
 852         gvdwc(l,k)=gvdwc(l,k)+gg(l)
 853       enddo
 854 C now derivative over rjk
 855       fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
 856      &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
 857             gg(1)=xjk*fac/rjk
 858             gg(2)=yjk*fac/rjk
 859             gg(3)=zjk*fac/rjk
 860       do m=1,3
 861         gvdwx(m,j)=gvdwx(m,j)-gg(m)
 862         gvdwx(m,k)=gvdwx(m,k)+gg(m)
 863       enddo
 864       do l=1,3
 865         gvdwc(l,j)=gvdwc(l,j)-gg(l)
 866         gvdwc(l,k)=gvdwc(l,k)+gg(l)
 867       enddo
 868       return
 869       end