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       nssbond=0
 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 (iout,'(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 c      write (iout,'(2(a,i5),4(a,f7.2))') "resi",resi," resj",resj,
 280 c     &  " ljxm",ljxm," ljxs",ljxs," ssxm",ssxm," rij",rij
 281       if (rij.gt.ljxm) then
 282         havebond=.false.
 283         ljd=rij-ljXs
 284         fac=(1.0D0/ljd)**expon
 285         e1=fac*fac*aa
 286         e2=fac*bb
 287         eij=eps1*eps2rt*eps3rt*(e1+e2)
 288         eps2der=eij*eps3rt
 289         eps3der=eij*eps2rt
 290         eij=eij*eps2rt*eps3rt*sss
 291
 292         sigder=-sig/sigsq
 293         e1=e1*eps1*eps2rt**2*eps3rt**2
 294         ed=-expon*(e1+eij)/ljd
 295         sigder=ed*sigder
 296         ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
 297         eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
 298         eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
 299         eom12=eij*eps1_om12+eps2der*eps2rt_om12
 300      &       -2.0D0*alf12*eps3der+sigder*sigsq_om12
 301       else if (rij.lt.ssxm) then
 302         havebond=.true.
 303         nssbond=nssbond+1
 304 c        write (iout,*) "ssMD: nssbond",nssbond
 305         ssd=rij-ssXs
 306         eij=ssA*ssd*ssd+ssB*ssd+ssC
 307         eij=eij*sss
 308         ed=2*akcm*ssd+akct*deltat12
 309         ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
 310         pom1=akct*ssd
 311         pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
 312         eom1=-2*akth*deltat1-pom1-om2*pom2
 313         eom2= 2*akth*deltat2+pom1-om1*pom2
 314         eom12=pom2
 315       else
 316 c          nssbond=nssbond+1
 317         omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
 318
 319         d_ssxm(1)=0.5D0*akct/ssA
 320         d_ssxm(2)=-d_ssxm(1)
 321         d_ssxm(3)=0.0D0
 322
 323         d_ljxm(1)=sig0ij/sqrt(sigsq**3)
 324         d_ljxm(2)=d_ljxm(1)*sigsq_om2
 325         d_ljxm(3)=d_ljxm(1)*sigsq_om12
 326         d_ljxm(1)=d_ljxm(1)*sigsq_om1
 327
 328 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
 329         xm=0.5d0*(ssxm+ljxm)
 330         do k=1,3
 331           d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
 332         enddo
 333         if (rij.lt.xm) then
 334           havebond=.true.
 335           ssm=ssC-0.25D0*ssB*ssB/ssA
 336           d_ssm(1)=0.5D0*akct*ssB/ssA
 337           d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
 338           d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
 339           d_ssm(3)=omega
 340           f1=(rij-xm)/(ssxm-xm)
 341           f2=(rij-ssxm)/(xm-ssxm)
 342           h1=h_base(f1,hd1)
 343           h2=h_base(f2,hd2)
 344           eij=ssm*h1+Ht*h2
 345           delta_inv=1.0d0/(xm-ssxm)
 346           deltasq_inv=delta_inv*delta_inv
 347           fac=ssm*hd1-Ht*hd2
 348           fac1=deltasq_inv*fac*(xm-rij)
 349           fac2=deltasq_inv*fac*(rij-ssxm)
 350           ed=delta_inv*(Ht*hd2-ssm*hd1)
 351           eij=eij*sss
 352           ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
 353           eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
 354           eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
 355           eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
 356         else
 357           havebond=.false.
 358           ljm=-0.25D0*ljB*bb/aa
 359           d_ljm(1)=-0.5D0*bb/aa*ljB
 360           d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
 361           d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
 362      +         alf12/eps3rt)
 363           d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
 364           f1=(rij-ljxm)/(xm-ljxm)
 365           f2=(rij-xm)/(ljxm-xm)
 366           h1=h_base(f1,hd1)
 367           h2=h_base(f2,hd2)
 368           eij=Ht*h1+ljm*h2
 369           delta_inv=1.0d0/(ljxm-xm)
 370           deltasq_inv=delta_inv*delta_inv
 371           fac=Ht*hd1-ljm*hd2
 372           fac1=deltasq_inv*fac*(ljxm-rij)
 373           fac2=deltasq_inv*fac*(rij-xm)
 374           ed=delta_inv*(ljm*hd2-Ht*hd1)
 375           eij=eij*sss
 376           ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
 377           eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
 378           eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
 379           eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
 380         endif
 381 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
 382
 383 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 384 c$$$        ssd=rij-ssXs
 385 c$$$        ljd=rij-ljXs
 386 c$$$        fac1=rij-ljxm
 387 c$$$        fac2=rij-ssxm
 388 c$$$
 389 c$$$        d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
 390 c$$$        d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
 391 c$$$        d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
 392 c$$$
 393 c$$$        ssm=ssC-0.25D0*ssB*ssB/ssA
 394 c$$$        d_ssm(1)=0.5D0*akct*ssB/ssA
 395 c$$$        d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
 396 c$$$        d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
 397 c$$$        d_ssm(3)=omega
 398 c$$$
 399 c$$$        ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
 400 c$$$        do k=1,3
 401 c$$$          d_ljm(k)=ljm*d_ljB(k)
 402 c$$$        enddo
 403 c$$$        ljm=ljm*ljB
 404 c$$$
 405 c$$$        ss=ssA*ssd*ssd+ssB*ssd+ssC
 406 c$$$        d_ss(0)=2.0d0*ssA*ssd+ssB
 407 c$$$        d_ss(2)=akct*ssd
 408 c$$$        d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
 409 c$$$        d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
 410 c$$$        d_ss(3)=omega
 411 c$$$
 412 c$$$        ljf=bb(itypi,itypj)/aa(itypi,itypj)
 413 c$$$        ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
 414 c$$$        d_ljf(0)=ljf*2.0d0*ljB*fac1
 415 c$$$        do k=1,3
 416 c$$$          d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
 417 c$$$     &         2.0d0*ljB*fac1*d_ljxm(k))
 418 c$$$        enddo
 419 c$$$        ljf=ljm+ljf*ljB*fac1*fac1
 420 c$$$
 421 c$$$        f1=(rij-ljxm)/(ssxm-ljxm)
 422 c$$$        f2=(rij-ssxm)/(ljxm-ssxm)
 423 c$$$        h1=h_base(f1,hd1)
 424 c$$$        h2=h_base(f2,hd2)
 425 c$$$        eij=ss*h1+ljf*h2
 426 c$$$        delta_inv=1.0d0/(ljxm-ssxm)
 427 c$$$        deltasq_inv=delta_inv*delta_inv
 428 c$$$        fac=ljf*hd2-ss*hd1
 429 c$$$        ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
 430 c$$$        eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
 431 c$$$     &       (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
 432 c$$$        eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
 433 c$$$     &       (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
 434 c$$$        eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
 435 c$$$     &       (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
 436 c$$$
 437 c$$$        havebond=.false.
 438 c$$$        if (ed.gt.0.0d0) havebond=.true.
 439 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 440
 441       endif
 442
 443       if (havebond) then
 444 #ifndef CLUST
 445 #ifndef WHAM
 446 c        if (dyn_ssbond_ij(i,j).eq.1.0d300) then
 447 c          write(iout,'(a15,f12.2,f8.1,2i5)')
 448 c     &         "SSBOND_E_FORM",totT,t_bath,i,j
 449 c        endif
 450 #endif
 451 #endif
 452         dyn_ssbond_ij(ici,icj)=eij
 453       else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300)
 454      &then
 455         dyn_ssbond_ij(ici,icj)=1.0d300
 456 #ifndef CLUST
 457 #ifndef WHAM
 458 c        write(iout,'(a15,f12.2,f8.1,2i5)')
 459 c     &       "SSBOND_E_BREAK",totT,t_bath,i,j
 460 #endif
 461 #endif
 462       endif
 463
 464 c-------TESTING CODE
 465       if (checkstop) then
 466         if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
 467      &       "CHECKSTOP",rij,eij,ed
 468         echeck(jcheck)=eij
 469       endif
 470       enddo
 471       if (checkstop) then
 472         write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
 473       endif
 474       enddo
 475       if (checkstop) then
 476         transgrad=.true.
 477         checkstop=.false.
 478       endif
 479 c-------END TESTING CODE
 480             gg_lipi(3)=ssgradlipi*eij
 481             gg_lipj(3)=ssgradlipj*eij
 482
 483       do k=1,3
 484         dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
 485         dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
 486       enddo
 487       do k=1,3
 488         gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
 489       enddo
 490       do k=1,3
 491         gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
 492      &       +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
 493      &       +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 494         gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
 495      &       +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
 496      &       +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
 497       enddo
 498 cgrad      do k=i,j-1
 499 cgrad        do l=1,3
 500 cgrad          gvdwc(l,k)=gvdwc(l,k)+gg(l)
 501 cgrad        enddo
 502 cgrad      enddo
 503
 504       do l=1,3
 505         gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
 506         gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
 507       enddo
 508
 509       return
 510       end
 511 C-----------------------------------------------------------------------------
 512
 513       double precision function h_base(x,deriv)
 514 c     A smooth function going 0->1 in range [0,1]
 515 c     It should NOT be called outside range [0,1], it will not work there.
 516       implicit none
 517
 518 c     Input arguments
 519       double precision x
 520
 521 c     Output arguments
 522       double precision deriv
 523
 524 c     Local variables
 525       double precision xsq
 526
 527
 528 c     Two parabolas put together.  First derivative zero at extrema
 529 c$$$      if (x.lt.0.5D0) then
 530 c$$$        h_base=2.0D0*x*x
 531 c$$$        deriv=4.0D0*x
 532 c$$$      else
 533 c$$$        deriv=1.0D0-x
 534 c$$$        h_base=1.0D0-2.0D0*deriv*deriv
 535 c$$$        deriv=4.0D0*deriv
 536 c$$$      endif
 537
 538 c     Third degree polynomial.  First derivative zero at extrema
 539       h_base=x*x*(3.0d0-2.0d0*x)
 540       deriv=6.0d0*x*(1.0d0-x)
 541
 542 c     Fifth degree polynomial.  First and second derivatives zero at extrema
 543 c$$$      xsq=x*x
 544 c$$$      h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
 545 c$$$      deriv=x-1.0d0
 546 c$$$      deriv=deriv*deriv
 547 c$$$      deriv=30.0d0*xsq*deriv
 548
 549       return
 550       end
 551 c----------------------------------------------------------------------------
 552       subroutine dyn_set_nss
 553 c     Adjust nss and other relevant variables based on dyn_ssbond_ij
 554 c      implicit none
 555
 556 c     Includes
 557       include 'DIMENSIONS'
 558 #ifdef MPI
 559       include "mpif.h"
 560 #endif
 561       include 'COMMON.SBRIDGE'
 562       include 'COMMON.CHAIN'
 563       include 'COMMON.IOUNITS'
 564 #ifndef CLUST
 565 #ifndef WHAM
 566 C      include 'COMMON.MD'
 567 #endif
 568 #endif
 569
 570 c     Local variables
 571       double precision emin
 572       integer i,j,imin
 573       integer diff,allflag(maxdim_cont),allnss,
 574      &     allihpb(maxdim_cont),alljhpb(maxdim_cont),
 575      &     newnss,newihpb(maxdim_cont),newjhpb(maxdim_cont)
 576       logical found
 577       integer i_newnss(1024),displ(0:1024)
 578       integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss
 579       nfgtasks=1
 580
 581       allnss=0
 582       do i=1,ns-1
 583         do j=i+1,ns
 584           if (dyn_ssbond_ij(i,j).lt.1.0d300) then
 585             allnss=allnss+1
 586             allflag(allnss)=0
 587             allihpb(allnss)=i
 588             alljhpb(allnss)=j
 589           endif
 590         enddo
 591       enddo
 592
 593 cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
 594
 595  1    emin=1.0d300
 596       do i=1,allnss
 597         if (allflag(i).eq.0 .and.
 598      &       dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
 599           emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
 600           imin=i
 601         endif
 602       enddo
 603       if (emin.lt.1.0d300) then
 604         allflag(imin)=1
 605         do i=1,allnss
 606           if (allflag(i).eq.0 .and.
 607      &         (allihpb(i).eq.allihpb(imin) .or.
 608      &         alljhpb(i).eq.allihpb(imin) .or.
 609      &         allihpb(i).eq.alljhpb(imin) .or.
 610      &         alljhpb(i).eq.alljhpb(imin))) then
 611             allflag(i)=-1
 612           endif
 613         enddo
 614         goto 1
 615       endif
 616
 617 cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
 618
 619       newnss=0
 620       do i=1,allnss
 621         if (allflag(i).eq.1) then
 622           newnss=newnss+1
 623           newihpb(newnss)=allihpb(i)
 624           newjhpb(newnss)=alljhpb(i)
 625         endif
 626       enddo
 627
 628 #ifdef MPI
 629       if (nfgtasks.gt.1)then
 630
 631         call MPI_Reduce(newnss,g_newnss,1,
 632      &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
 633         call MPI_Gather(newnss,1,MPI_INTEGER,
 634      &                  i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
 635         displ(0)=0
 636         do i=1,nfgtasks-1,1
 637           displ(i)=i_newnss(i-1)+displ(i-1)
 638         enddo
 639         call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
 640      &                   g_newihpb,i_newnss,displ,MPI_INTEGER,
 641      &                   king,FG_COMM,IERR)
 642         call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
 643      &                   g_newjhpb,i_newnss,displ,MPI_INTEGER,
 644      &                   king,FG_COMM,IERR)
 645         if(fg_rank.eq.0) then
 646 c         print *,'g_newnss',g_newnss
 647 c         print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
 648 c         print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
 649          newnss=g_newnss
 650          do i=1,newnss
 651           newihpb(i)=g_newihpb(i)
 652           newjhpb(i)=g_newjhpb(i)
 653          enddo
 654         endif
 655       endif
 656 #endif
 657
 658       diff=newnss-nss
 659
 660 cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
 661
 662       do i=1,nss
 663         found=.false.
 664         do j=1,newnss
 665           if (idssb(i).eq.newihpb(j) .and.
 666      &         jdssb(i).eq.newjhpb(j)) found=.true.
 667         enddo
 668 #ifndef CLUST
 669 #ifndef WHAM
 670 c        if (.not.found.and.fg_rank.eq.0)
 671 c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
 672 c     &       "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
 673 #endif
 674 #endif
 675       enddo
 676
 677       do i=1,newnss
 678         found=.false.
 679         do j=1,nss
 680           if (newihpb(i).eq.idssb(j) .and.
 681      &         newjhpb(i).eq.jdssb(j)) found=.true.
 682         enddo
 683 #ifndef CLUST
 684 #ifndef WHAM
 685 c        if (.not.found.and.fg_rank.eq.0)
 686 c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
 687 c     &       "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
 688 #endif
 689 #endif
 690       enddo
 691
 692       nss=newnss
 693       do i=1,nss
 694         idssb(i)=newihpb(i)
 695         jdssb(i)=newjhpb(i)
 696       enddo
 697
 698       return
 699       end
 700
 701 c----------------------------------------------------------------------------
 702
 703 #ifdef SSREAD
 704 #ifdef WHAM
 705       subroutine read_ssHist
 706       implicit none
 707
 708 c     Includes
 709       include 'DIMENSIONS'
 710       include "DIMENSIONS.FREE"
 711       include 'COMMON.FREE'
 712
 713 c     Local variables
 714       integer i,j
 715       character*80 controlcard
 716
 717       do i=1,dyn_nssHist
 718         call card_concat(controlcard,.true.)
 719         read(controlcard,*)
 720      &       dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
 721       enddo
 722
 723       return
 724       end
 725 #endif
 726 #endif
 727 c$$$C----------------------------------------------------------------------------
 728       subroutine triple_ssbond_ene(resi,resj,resk,eij)
 729       include 'DIMENSIONS'
 730       include 'COMMON.SBRIDGE'
 731       include 'COMMON.CHAIN'
 732       include 'COMMON.DERIV'
 733       include 'COMMON.LOCAL'
 734       include 'COMMON.INTERACT'
 735       include 'COMMON.VAR'
 736       include 'COMMON.IOUNITS'
 737       include 'COMMON.CALC'
 738 #ifndef CLUST
 739 #ifndef WHAM
 740 C      include 'COMMON.MD'
 741 #endif
 742 #endif
 743
 744 c     External functions
 745       double precision h_base
 746       external h_base
 747
 748 c     Input arguments
 749       integer resi,resj,resk
 750
 751 c     Output arguments
 752       double precision eij,eij1,eij2,eij3
 753
 754 c     Local variables
 755       logical havebond
 756 c      integer itypi,itypj,k,l
 757       double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
 758       double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
 759       double precision xik,yik,zik,xjk,yjk,zjk
 760       double precision sig0ij,ljd,sig,fac,e1,e2
 761       double precision dcosom1(3),dcosom2(3),ed
 762       double precision pom1,pom2
 763       double precision ljA,ljB,ljXs
 764       double precision d_ljB(1:3)
 765       double precision ssA,ssB,ssC,ssXs
 766       double precision ssxm,ljxm,ssm,ljm
 767       double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
 768
 769       i=resi
 770       j=resj
 771       k=resk
 772 C      write(iout,*) resi,resj,resk
 773       itypi=itype(i)
 774       dxi=dc_norm(1,nres+i)
 775       dyi=dc_norm(2,nres+i)
 776       dzi=dc_norm(3,nres+i)
 777       dsci_inv=vbld_inv(i+nres)
 778       xi=c(1,nres+i)
 779       yi=c(2,nres+i)
 780       zi=c(3,nres+i)
 781
 782       itypj=itype(j)
 783       xj=c(1,nres+j)
 784       yj=c(2,nres+j)
 785       zj=c(3,nres+j)
 786
 787       dxj=dc_norm(1,nres+j)
 788       dyj=dc_norm(2,nres+j)
 789       dzj=dc_norm(3,nres+j)
 790       dscj_inv=vbld_inv(j+nres)
 791       itypk=itype(k)
 792       xk=c(1,nres+k)
 793       yk=c(2,nres+k)
 794       zk=c(3,nres+k)
 795
 796       dxk=dc_norm(1,nres+k)
 797       dyk=dc_norm(2,nres+k)
 798       dzk=dc_norm(3,nres+k)
 799       dscj_inv=vbld_inv(k+nres)
 800       xij=xj-xi
 801       xik=xk-xi
 802       xjk=xk-xj
 803       yij=yj-yi
 804       yik=yk-yi
 805       yjk=yk-yj
 806       zij=zj-zi
 807       zik=zk-zi
 808       zjk=zk-zj
 809       rrij=(xij*xij+yij*yij+zij*zij)
 810       rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
 811       rrik=(xik*xik+yik*yik+zik*zik)
 812       rik=dsqrt(rrik)
 813       rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
 814       rjk=dsqrt(rrjk)
 815 C there are three combination of distances for each trisulfide bonds
 816 C The first case the ith atom is the center
 817 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
 818 C distance y is second distance the a,b,c,d are parameters derived for
 819 C this problem d parameter was set as a penalty currenlty set to 1.
 820       eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
 821 C second case jth atom is center
 822       eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
 823 C the third case kth atom is the center
 824       eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
 825 C      eij2=0.0
 826 C      eij3=0.0
 827 C      eij1=0.0
 828       eij=eij1+eij2+eij3
 829 C      write(iout,*)i,j,k,eij
 830 C The energy penalty calculated now time for the gradient part
 831 C derivative over rij
 832       fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
 833      &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
 834             gg(1)=xij*fac/rij
 835             gg(2)=yij*fac/rij
 836             gg(3)=zij*fac/rij
 837       do m=1,3
 838         gvdwx(m,i)=gvdwx(m,i)-gg(m)
 839         gvdwx(m,j)=gvdwx(m,j)+gg(m)
 840       enddo
 841       do l=1,3
 842         gvdwc(l,i)=gvdwc(l,i)-gg(l)
 843         gvdwc(l,j)=gvdwc(l,j)+gg(l)
 844       enddo
 845 C now derivative over rik
 846       fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
 847      &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
 848             gg(1)=xik*fac/rik
 849             gg(2)=yik*fac/rik
 850             gg(3)=zik*fac/rik
 851       do m=1,3
 852         gvdwx(m,i)=gvdwx(m,i)-gg(m)
 853         gvdwx(m,k)=gvdwx(m,k)+gg(m)
 854       enddo
 855       do l=1,3
 856         gvdwc(l,i)=gvdwc(l,i)-gg(l)
 857         gvdwc(l,k)=gvdwc(l,k)+gg(l)
 858       enddo
 859 C now derivative over rjk
 860       fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
 861      &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
 862             gg(1)=xjk*fac/rjk
 863             gg(2)=yjk*fac/rjk
 864             gg(3)=zjk*fac/rjk
 865       do m=1,3
 866         gvdwx(m,j)=gvdwx(m,j)-gg(m)
 867         gvdwx(m,k)=gvdwx(m,k)+gg(m)
 868       enddo
 869       do l=1,3
 870         gvdwc(l,j)=gvdwc(l,j)-gg(l)
 871         gvdwc(l,k)=gvdwc(l,k)+gg(l)
 872       enddo
 873       return
 874       end