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