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