-c----------------------------------------------------------------------------
-
-
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-
-c$$$c-----------------------------------------------------------------------------
-c$$$
-c$$$ subroutine ss_relax(i_in,j_in)
-c$$$ implicit none
-c$$$
-c$$$c Includes
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.CHAIN'
-c$$$ include 'COMMON.IOUNITS'
-c$$$ include 'COMMON.INTERACT'
-c$$$
-c$$$c Input arguments
-c$$$ integer i_in,j_in
-c$$$
-c$$$c Local variables
-c$$$ integer i,iretcode,nfun_sc
-c$$$ logical scfail
-c$$$ double precision var(maxvar),e_sc,etot
-c$$$
-c$$$
-c$$$ mask_r=.true.
-c$$$ do i=nnt,nct
-c$$$ mask_side(i)=0
-c$$$ enddo
-c$$$ mask_side(i_in)=1
-c$$$ mask_side(j_in)=1
-c$$$
-c$$$c Minimize the two selected side-chains
-c$$$ call overlap_sc(scfail) ! Better not fail!
-c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
-c$$$
-c$$$ mask_r=.false.
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-c$$$c-------------------------------------------------------------
-c$$$
-c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
-c$$$c Minimize side-chains only, starting from geom but without modifying
-c$$$c bond lengths.
-c$$$c If mask_r is already set, only the selected side-chains are minimized,
-c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
-c$$$ implicit none
-c$$$
-c$$$c Includes
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.IOUNITS'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.CHAIN'
-c$$$ include 'COMMON.GEO'
-c$$$ include 'COMMON.MINIM'
-c$$$ integer icall
-c$$$ common /srutu/ icall
-c$$$
-c$$$c Output arguments
-c$$$ double precision etot_sc
-c$$$ integer iretcode,nfun
-c$$$
-c$$$c External functions/subroutines
-c$$$ external func_sc,grad_sc,fdum
-c$$$
-c$$$c Local variables
-c$$$ integer liv,lv
-c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
-c$$$ integer iv(liv)
-c$$$ double precision rdum(1)
-c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
-c$$$ integer idum(1)
-c$$$ integer i,nvar_restr
-c$$$
-c$$$
-c$$$cmc start_minim=.true.
-c$$$ call deflt(2,iv,liv,lv,v)
-c$$$* 12 means fresh start, dont call deflt
-c$$$ iv(1)=12
-c$$$* max num of fun calls
-c$$$ if (maxfun.eq.0) maxfun=500
-c$$$ iv(17)=maxfun
-c$$$* max num of iterations
-c$$$ if (maxmin.eq.0) maxmin=1000
-c$$$ iv(18)=maxmin
-c$$$* controls output
-c$$$ iv(19)=1
-c$$$* selects output unit
-c$$$ iv(21)=0
-c$$$c iv(21)=iout ! DEBUG
-c$$$c iv(21)=8 ! DEBUG
-c$$$* 1 means to print out result
-c$$$ iv(22)=0
-c$$$c iv(22)=1 ! DEBUG
-c$$$* 1 means to print out summary stats
-c$$$ iv(23)=0
-c$$$c iv(23)=1 ! DEBUG
-c$$$* 1 means to print initial x and d
-c$$$ iv(24)=0
-c$$$c iv(24)=1 ! DEBUG
-c$$$* min val for v(radfac) default is 0.1
-c$$$ v(24)=0.1D0
-c$$$* max val for v(radfac) default is 4.0
-c$$$ v(25)=2.0D0
-c$$$c v(25)=4.0D0
-c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
-c$$$* the sumsl default is 0.1
-c$$$ v(26)=0.1D0
-c$$$* false conv if (act fnctn decrease) .lt. v(34)
-c$$$* the sumsl default is 100*machep
-c$$$ v(34)=v(34)/100.0D0
-c$$$* absolute convergence
-c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
-c$$$ v(31)=tolf
-c$$$* relative convergence
-c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
-c$$$ v(32)=rtolf
-c$$$* controls initial step size
-c$$$ v(35)=1.0D-1
-c$$$* large vals of d correspond to small components of step
-c$$$ do i=1,nphi
-c$$$ d(i)=1.0D-1
-c$$$ enddo
-c$$$ do i=nphi+1,nvar
-c$$$ d(i)=1.0D-1
-c$$$ enddo
-c$$$
-c$$$ call geom_to_var(nvar,x)
-c$$$ IF (mask_r) THEN
-c$$$ do i=1,nres ! Just in case...
-c$$$ mask_phi(i)=0
-c$$$ mask_theta(i)=0
-c$$$ enddo
-c$$$ call x2xx(x,xx,nvar_restr)
-c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
-c$$$ & iv,liv,lv,v,idum,rdum,fdum)
-c$$$ call xx2x(x,xx)
-c$$$ ELSE
-c$$$c When minimizing ALL side-chains, etotal_sc is a little
-c$$$c faster if we don't set mask_r
-c$$$ do i=1,nres
-c$$$ mask_phi(i)=0
-c$$$ mask_theta(i)=0
-c$$$ mask_side(i)=1
-c$$$ enddo
-c$$$ call x2xx(x,xx,nvar_restr)
-c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
-c$$$ & iv,liv,lv,v,idum,rdum,fdum)
-c$$$ call xx2x(x,xx)
-c$$$ ENDIF
-c$$$ call var_to_geom(nvar,x)
-c$$$ call chainbuild_sc
-c$$$ etot_sc=v(10)
-c$$$ iretcode=iv(1)
-c$$$ nfun=iv(6)
-c$$$ return
-c$$$ end
-c$$$
-c$$$C--------------------------------------------------------------------------
-c$$$
-c$$$ subroutine chainbuild_sc
-c$$$ implicit none
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.INTERACT'
-c$$$
-c$$$c Local variables
-c$$$ integer i
-c$$$
-c$$$
-c$$$ do i=nnt,nct
-c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
-c$$$ call locate_side_chain(i)
-c$$$ endif
-c$$$ enddo
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-c$$$C--------------------------------------------------------------------------
-c$$$
-c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
-c$$$ implicit none
-c$$$
-c$$$c Includes
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.DERIV'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.MINIM'
-c$$$ include 'COMMON.IOUNITS'
-c$$$
-c$$$c Input arguments
-c$$$ integer n
-c$$$ double precision x(maxvar)
-c$$$ double precision ufparm
-c$$$ external ufparm
-c$$$
-c$$$c Input/Output arguments
-c$$$ integer nf
-c$$$ integer uiparm(1)
-c$$$ double precision urparm(1)
-c$$$
-c$$$c Output arguments
-c$$$ double precision f
-c$$$
-c$$$c Local variables
-c$$$ double precision energia(0:n_ene)
-c$$$#ifdef OSF
-c$$$c Variables used to intercept NaNs
-c$$$ double precision x_sum
-c$$$ integer i_NAN
-c$$$#endif
-c$$$
-c$$$
-c$$$ nfl=nf
-c$$$ icg=mod(nf,2)+1
-c$$$
-c$$$#ifdef OSF
-c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
-c$$$ x_sum=0.D0
-c$$$ do i_NAN=1,n
-c$$$ x_sum=x_sum+x(i_NAN)
-c$$$ enddo
-c$$$c Calculate the energy only if the coordinates are ok
-c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
-c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
-c$$$ f=1.0D+77
-c$$$ nf=0
-c$$$ else
-c$$$#endif
-c$$$
-c$$$ call var_to_geom_restr(n,x)
-c$$$ call zerograd
-c$$$ call chainbuild_sc
-c$$$ call etotal_sc(energia(0))
-c$$$ f=energia(0)
-c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
-c$$$
-c$$$#ifdef OSF
-c$$$ endif
-c$$$#endif
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-c$$$c-------------------------------------------------------
-c$$$
-c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
-c$$$ implicit none
-c$$$
-c$$$c Includes
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.CHAIN'
-c$$$ include 'COMMON.DERIV'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.INTERACT'
-c$$$ include 'COMMON.MINIM'
-c$$$
-c$$$c Input arguments
-c$$$ integer n
-c$$$ double precision x(maxvar)
-c$$$ double precision ufparm
-c$$$ external ufparm
-c$$$
-c$$$c Input/Output arguments
-c$$$ integer nf
-c$$$ integer uiparm(1)
-c$$$ double precision urparm(1)
-c$$$
-c$$$c Output arguments
-c$$$ double precision g(maxvar)
-c$$$
-c$$$c Local variables
-c$$$ double precision f,gphii,gthetai,galphai,gomegai
-c$$$ integer ig,ind,i,j,k,igall,ij
-c$$$
-c$$$
-c$$$ icg=mod(nf,2)+1
-c$$$ if (nf-nfl+1) 20,30,40
-c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
-c$$$c write (iout,*) 'grad 20'
-c$$$ if (nf.eq.0) return
-c$$$ goto 40
-c$$$ 30 call var_to_geom_restr(n,x)
-c$$$ call chainbuild_sc
-c$$$C
-c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
-c$$$C
-c$$$ 40 call cartder
-c$$$C
-c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
-c$$$C
-c$$$
-c$$$ ig=0
-c$$$ ind=nres-2
-c$$$ do i=2,nres-2
-c$$$ IF (mask_phi(i+2).eq.1) THEN
-c$$$ gphii=0.0D0
-c$$$ do j=i+1,nres-1
-c$$$ ind=ind+1
-c$$$ do k=1,3
-c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
-c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
-c$$$ enddo
-c$$$ enddo
-c$$$ ig=ig+1
-c$$$ g(ig)=gphii
-c$$$ ELSE
-c$$$ ind=ind+nres-1-i
-c$$$ ENDIF
-c$$$ enddo
-c$$$
-c$$$
-c$$$ ind=0
-c$$$ do i=1,nres-2
-c$$$ IF (mask_theta(i+2).eq.1) THEN
-c$$$ ig=ig+1
-c$$$ gthetai=0.0D0
-c$$$ do j=i+1,nres-1
-c$$$ ind=ind+1
-c$$$ do k=1,3
-c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
-c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
-c$$$ enddo
-c$$$ enddo
-c$$$ g(ig)=gthetai
-c$$$ ELSE
-c$$$ ind=ind+nres-1-i
-c$$$ ENDIF
-c$$$ enddo
-c$$$
-c$$$ do i=2,nres-1
-c$$$ if (itype(i).ne.10) then
-c$$$ IF (mask_side(i).eq.1) THEN
-c$$$ ig=ig+1
-c$$$ galphai=0.0D0
-c$$$ do k=1,3
-c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
-c$$$ enddo
-c$$$ g(ig)=galphai
-c$$$ ENDIF
-c$$$ endif
-c$$$ enddo
-c$$$
-c$$$
-c$$$ do i=2,nres-1
-c$$$ if (itype(i).ne.10) then
-c$$$ IF (mask_side(i).eq.1) THEN
-c$$$ ig=ig+1
-c$$$ gomegai=0.0D0
-c$$$ do k=1,3
-c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
-c$$$ enddo
-c$$$ g(ig)=gomegai
-c$$$ ENDIF
-c$$$ endif
-c$$$ enddo
-c$$$
-c$$$C
-c$$$C Add the components corresponding to local energy terms.
-c$$$C
-c$$$
-c$$$ ig=0
-c$$$ igall=0
-c$$$ do i=4,nres
-c$$$ igall=igall+1
-c$$$ if (mask_phi(i).eq.1) then
-c$$$ ig=ig+1
-c$$$ g(ig)=g(ig)+gloc(igall,icg)
-c$$$ endif
-c$$$ enddo
-c$$$
-c$$$ do i=3,nres
-c$$$ igall=igall+1
-c$$$ if (mask_theta(i).eq.1) then
-c$$$ ig=ig+1
-c$$$ g(ig)=g(ig)+gloc(igall,icg)
-c$$$ endif
-c$$$ enddo
-c$$$
-c$$$ do ij=1,2
-c$$$ do i=2,nres-1
-c$$$ if (itype(i).ne.10) then
-c$$$ igall=igall+1
-c$$$ if (mask_side(i).eq.1) then
-c$$$ ig=ig+1
-c$$$ g(ig)=g(ig)+gloc(igall,icg)
-c$$$ endif
-c$$$ endif
-c$$$ enddo
-c$$$ enddo
-c$$$
-c$$$cd do i=1,ig
-c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
-c$$$cd enddo
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-c$$$C-----------------------------------------------------------------------------
-c$$$
-c$$$ subroutine etotal_sc(energy_sc)
-c$$$ implicit none
-c$$$
-c$$$c Includes
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.INTERACT'
-c$$$ include 'COMMON.DERIV'
-c$$$ include 'COMMON.FFIELD'
-c$$$
-c$$$c Output arguments
-c$$$ double precision energy_sc(0:n_ene)
-c$$$
-c$$$c Local variables
-c$$$ double precision evdw,escloc
-c$$$ integer i,j
-c$$$
-c$$$
-c$$$ do i=1,n_ene
-c$$$ energy_sc(i)=0.0D0
-c$$$ enddo
-c$$$
-c$$$ if (mask_r) then
-c$$$ call egb_sc(evdw)
-c$$$ call esc_sc(escloc)
-c$$$ else
-c$$$ call egb(evdw)
-c$$$ call esc(escloc)
-c$$$ endif
-c$$$
-c$$$ if (evdw.eq.1.0D20) then
-c$$$ energy_sc(0)=evdw
-c$$$ else
-c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
-c$$$ endif
-c$$$ energy_sc(1)=evdw
-c$$$ energy_sc(12)=escloc
-c$$$
-c$$$C
-c$$$C Sum up the components of the Cartesian gradient.
-c$$$C
-c$$$ do i=1,nct
-c$$$ do j=1,3
-c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
-c$$$ enddo
-c$$$ enddo
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-c$$$C-----------------------------------------------------------------------------
-c$$$
-c$$$ subroutine egb_sc(evdw)
-c$$$C
-c$$$C This subroutine calculates the interaction energy of nonbonded side chains
-c$$$C assuming the Gay-Berne potential of interaction.
-c$$$C
-c$$$ implicit real*8 (a-h,o-z)
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.GEO'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.LOCAL'
-c$$$ include 'COMMON.CHAIN'
-c$$$ include 'COMMON.DERIV'
-c$$$ include 'COMMON.NAMES'
-c$$$ include 'COMMON.INTERACT'
-c$$$ include 'COMMON.IOUNITS'
-c$$$ include 'COMMON.CALC'
-c$$$ include 'COMMON.CONTROL'
-c$$$ logical lprn
-c$$$ evdw=0.0D0
-c$$$ energy_dec=.false.
-c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
-c$$$ evdw=0.0D0
-c$$$ lprn=.false.
-c$$$c if (icall.eq.0) lprn=.false.
-c$$$ ind=0
-c$$$ do i=iatsc_s,iatsc_e
-c$$$ itypi=itype(i)
-c$$$ itypi1=itype(i+1)
-c$$$ xi=c(1,nres+i)
-c$$$ yi=c(2,nres+i)
-c$$$ zi=c(3,nres+i)
-c$$$ dxi=dc_norm(1,nres+i)
-c$$$ dyi=dc_norm(2,nres+i)
-c$$$ dzi=dc_norm(3,nres+i)
-c$$$c dsci_inv=dsc_inv(itypi)
-c$$$ dsci_inv=vbld_inv(i+nres)
-c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
-c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
-c$$$C
-c$$$C Calculate SC interaction energy.
-c$$$C
-c$$$ do iint=1,nint_gr(i)
-c$$$ do j=istart(i,iint),iend(i,iint)
-c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
-c$$$ ind=ind+1
-c$$$ itypj=itype(j)
-c$$$c dscj_inv=dsc_inv(itypj)
-c$$$ dscj_inv=vbld_inv(j+nres)
-c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
-c$$$c & 1.0d0/vbld(j+nres)
-c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
-c$$$ sig0ij=sigma(itypi,itypj)
-c$$$ chi1=chi(itypi,itypj)
-c$$$ chi2=chi(itypj,itypi)
-c$$$ chi12=chi1*chi2
-c$$$ chip1=chip(itypi)
-c$$$ chip2=chip(itypj)
-c$$$ chip12=chip1*chip2
-c$$$ alf1=alp(itypi)
-c$$$ alf2=alp(itypj)
-c$$$ alf12=0.5D0*(alf1+alf2)
-c$$$C For diagnostics only!!!
-c$$$c chi1=0.0D0
-c$$$c chi2=0.0D0
-c$$$c chi12=0.0D0
-c$$$c chip1=0.0D0
-c$$$c chip2=0.0D0
-c$$$c chip12=0.0D0
-c$$$c alf1=0.0D0
-c$$$c alf2=0.0D0
-c$$$c alf12=0.0D0
-c$$$ xj=c(1,nres+j)-xi
-c$$$ yj=c(2,nres+j)-yi
-c$$$ zj=c(3,nres+j)-zi
-c$$$ dxj=dc_norm(1,nres+j)
-c$$$ dyj=dc_norm(2,nres+j)
-c$$$ dzj=dc_norm(3,nres+j)
-c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
-c$$$c write (iout,*) "j",j," dc_norm",
-c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
-c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
-c$$$ rij=dsqrt(rrij)
-c$$$C Calculate angle-dependent terms of energy and contributions to their
-c$$$C derivatives.
-c$$$ call sc_angular
-c$$$ sigsq=1.0D0/sigsq
-c$$$ sig=sig0ij*dsqrt(sigsq)
-c$$$ rij_shift=1.0D0/rij-sig+sig0ij
-c$$$c for diagnostics; uncomment
-c$$$c rij_shift=1.2*sig0ij
-c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
-c$$$ if (rij_shift.le.0.0D0) then
-c$$$ evdw=1.0D20
-c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-c$$$cd & restyp(itypi),i,restyp(itypj),j,
-c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
-c$$$ return
-c$$$ endif
-c$$$ sigder=-sig*sigsq
-c$$$c---------------------------------------------------------------
-c$$$ rij_shift=1.0D0/rij_shift
-c$$$ fac=rij_shift**expon
-c$$$ e1=fac*fac*aa(itypi,itypj)
-c$$$ e2=fac*bb(itypi,itypj)
-c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
-c$$$ eps2der=evdwij*eps3rt
-c$$$ eps3der=evdwij*eps2rt
-c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
-c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
-c$$$ evdwij=evdwij*eps2rt*eps3rt
-c$$$ evdw=evdw+evdwij
-c$$$ if (lprn) then
-c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
-c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-c$$$ & restyp(itypi),i,restyp(itypj),j,
-c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
-c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
-c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
-c$$$ & evdwij
-c$$$ endif
-c$$$
-c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
-c$$$ & 'evdw',i,j,evdwij
-c$$$
-c$$$C Calculate gradient components.
-c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
-c$$$ fac=-expon*(e1+evdwij)*rij_shift
-c$$$ sigder=fac*sigder
-c$$$ fac=rij*fac
-c$$$c fac=0.0d0
-c$$$C Calculate the radial part of the gradient
-c$$$ gg(1)=xj*fac
-c$$$ gg(2)=yj*fac
-c$$$ gg(3)=zj*fac
-c$$$C Calculate angular part of the gradient.
-c$$$ call sc_grad
-c$$$ ENDIF
-c$$$ enddo ! j
-c$$$ enddo ! iint
-c$$$ enddo ! i
-c$$$ energy_dec=.false.
-c$$$ return
-c$$$ end
-c$$$
-c$$$c-----------------------------------------------------------------------------
-c$$$
-c$$$ subroutine esc_sc(escloc)
-c$$$C Calculate the local energy of a side chain and its derivatives in the
-c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
-c$$$C ALPHA and OMEGA.
-c$$$ implicit real*8 (a-h,o-z)
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.GEO'
-c$$$ include 'COMMON.LOCAL'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.INTERACT'
-c$$$ include 'COMMON.DERIV'
-c$$$ include 'COMMON.CHAIN'
-c$$$ include 'COMMON.IOUNITS'
-c$$$ include 'COMMON.NAMES'
-c$$$ include 'COMMON.FFIELD'
-c$$$ include 'COMMON.CONTROL'
-c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
-c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
-c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
-c$$$ delta=0.02d0*pi
-c$$$ escloc=0.0D0
-c$$$c write (iout,'(a)') 'ESC'
-c$$$ do i=loc_start,loc_end
-c$$$ IF (mask_side(i).eq.1) THEN
-c$$$ it=itype(i)
-c$$$ if (it.eq.10) goto 1
-c$$$ nlobit=nlob(it)
-c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
-c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
-c$$$ theti=theta(i+1)-pipol
-c$$$ x(1)=dtan(theti)
-c$$$ x(2)=alph(i)
-c$$$ x(3)=omeg(i)
-c$$$
-c$$$ if (x(2).gt.pi-delta) then
-c$$$ xtemp(1)=x(1)
-c$$$ xtemp(2)=pi-delta
-c$$$ xtemp(3)=x(3)
-c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
-c$$$ xtemp(2)=pi
-c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
-c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
-c$$$ & escloci,dersc(2))
-c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
-c$$$ & ddersc0(1),dersc(1))
-c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
-c$$$ & ddersc0(3),dersc(3))
-c$$$ xtemp(2)=pi-delta
-c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
-c$$$ xtemp(2)=pi
-c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
-c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
-c$$$ & dersc0(2),esclocbi,dersc02)
-c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
-c$$$ & dersc12,dersc01)
-c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
-c$$$ dersc0(1)=dersc01
-c$$$ dersc0(2)=dersc02
-c$$$ dersc0(3)=0.0d0
-c$$$ do k=1,3
-c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
-c$$$ enddo
-c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
-c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c$$$c & esclocbi,ss,ssd
-c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
-c$$$c escloci=esclocbi
-c$$$c write (iout,*) escloci
-c$$$ else if (x(2).lt.delta) then
-c$$$ xtemp(1)=x(1)
-c$$$ xtemp(2)=delta
-c$$$ xtemp(3)=x(3)
-c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
-c$$$ xtemp(2)=0.0d0
-c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
-c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
-c$$$ & escloci,dersc(2))
-c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
-c$$$ & ddersc0(1),dersc(1))
-c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
-c$$$ & ddersc0(3),dersc(3))
-c$$$ xtemp(2)=delta
-c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
-c$$$ xtemp(2)=0.0d0
-c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
-c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
-c$$$ & dersc0(2),esclocbi,dersc02)
-c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
-c$$$ & dersc12,dersc01)
-c$$$ dersc0(1)=dersc01
-c$$$ dersc0(2)=dersc02
-c$$$ dersc0(3)=0.0d0
-c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
-c$$$ do k=1,3
-c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
-c$$$ enddo
-c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
-c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c$$$c & esclocbi,ss,ssd
-c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
-c$$$c write (iout,*) escloci
-c$$$ else
-c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
-c$$$ endif
-c$$$
-c$$$ escloc=escloc+escloci
-c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
-c$$$ & 'escloc',i,escloci
-c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
-c$$$
-c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
-c$$$ & wscloc*dersc(1)
-c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
-c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
-c$$$ 1 continue
-c$$$ ENDIF
-c$$$ enddo
-c$$$ return
-c$$$ end
-c$$$
-c$$$C-----------------------------------------------------------------------------
-c$$$
-c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
-c$$$C
-c$$$C This subroutine calculates the interaction energy of nonbonded side chains
-c$$$C assuming the Gay-Berne potential of interaction.
-c$$$C
-c$$$ implicit real*8 (a-h,o-z)
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.GEO'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.LOCAL'
-c$$$ include 'COMMON.CHAIN'
-c$$$ include 'COMMON.DERIV'
-c$$$ include 'COMMON.NAMES'
-c$$$ include 'COMMON.INTERACT'
-c$$$ include 'COMMON.IOUNITS'
-c$$$ include 'COMMON.CALC'
-c$$$ include 'COMMON.CONTROL'
-c$$$ logical lprn
-c$$$ evdw=0.0D0
-c$$$ energy_dec=.false.
-c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
-c$$$ evdw=0.0D0
-c$$$ lprn=.false.
-c$$$ ind=0
-c$$$c$$$ do i=iatsc_s,iatsc_e
-c$$$ i=i_sc
-c$$$ itypi=itype(i)
-c$$$ itypi1=itype(i+1)
-c$$$ xi=c(1,nres+i)
-c$$$ yi=c(2,nres+i)
-c$$$ zi=c(3,nres+i)
-c$$$ dxi=dc_norm(1,nres+i)
-c$$$ dyi=dc_norm(2,nres+i)
-c$$$ dzi=dc_norm(3,nres+i)
-c$$$c dsci_inv=dsc_inv(itypi)
-c$$$ dsci_inv=vbld_inv(i+nres)
-c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
-c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
-c$$$C
-c$$$C Calculate SC interaction energy.
-c$$$C
-c$$$c$$$ do iint=1,nint_gr(i)
-c$$$c$$$ do j=istart(i,iint),iend(i,iint)
-c$$$ j=j_sc
-c$$$ ind=ind+1
-c$$$ itypj=itype(j)
-c$$$c dscj_inv=dsc_inv(itypj)
-c$$$ dscj_inv=vbld_inv(j+nres)
-c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
-c$$$c & 1.0d0/vbld(j+nres)
-c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
-c$$$ sig0ij=sigma(itypi,itypj)
-c$$$ chi1=chi(itypi,itypj)
-c$$$ chi2=chi(itypj,itypi)
-c$$$ chi12=chi1*chi2
-c$$$ chip1=chip(itypi)
-c$$$ chip2=chip(itypj)
-c$$$ chip12=chip1*chip2
-c$$$ alf1=alp(itypi)
-c$$$ alf2=alp(itypj)
-c$$$ alf12=0.5D0*(alf1+alf2)
-c$$$C For diagnostics only!!!
-c$$$c chi1=0.0D0
-c$$$c chi2=0.0D0
-c$$$c chi12=0.0D0
-c$$$c chip1=0.0D0
-c$$$c chip2=0.0D0
-c$$$c chip12=0.0D0
-c$$$c alf1=0.0D0
-c$$$c alf2=0.0D0
-c$$$c alf12=0.0D0
-c$$$ xj=c(1,nres+j)-xi
-c$$$ yj=c(2,nres+j)-yi
-c$$$ zj=c(3,nres+j)-zi
-c$$$ dxj=dc_norm(1,nres+j)
-c$$$ dyj=dc_norm(2,nres+j)
-c$$$ dzj=dc_norm(3,nres+j)
-c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
-c$$$c write (iout,*) "j",j," dc_norm",
-c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
-c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
-c$$$ rij=dsqrt(rrij)
-c$$$C Calculate angle-dependent terms of energy and contributions to their
-c$$$C derivatives.
-c$$$ call sc_angular
-c$$$ sigsq=1.0D0/sigsq
-c$$$ sig=sig0ij*dsqrt(sigsq)
-c$$$ rij_shift=1.0D0/rij-sig+sig0ij
-c$$$c for diagnostics; uncomment
-c$$$c rij_shift=1.2*sig0ij
-c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
-c$$$ if (rij_shift.le.0.0D0) then
-c$$$ evdw=1.0D20
-c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-c$$$cd & restyp(itypi),i,restyp(itypj),j,
-c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
-c$$$ return
-c$$$ endif
-c$$$ sigder=-sig*sigsq
-c$$$c---------------------------------------------------------------
-c$$$ rij_shift=1.0D0/rij_shift
-c$$$ fac=rij_shift**expon
-c$$$ e1=fac*fac*aa(itypi,itypj)
-c$$$ e2=fac*bb(itypi,itypj)
-c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
-c$$$ eps2der=evdwij*eps3rt
-c$$$ eps3der=evdwij*eps2rt
-c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
-c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
-c$$$ evdwij=evdwij*eps2rt*eps3rt
-c$$$ evdw=evdw+evdwij
-c$$$ if (lprn) then
-c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
-c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-c$$$ & restyp(itypi),i,restyp(itypj),j,
-c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
-c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
-c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
-c$$$ & evdwij
-c$$$ endif
-c$$$
-c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
-c$$$ & 'evdw',i,j,evdwij
-c$$$
-c$$$C Calculate gradient components.
-c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
-c$$$ fac=-expon*(e1+evdwij)*rij_shift
-c$$$ sigder=fac*sigder
-c$$$ fac=rij*fac
-c$$$c fac=0.0d0
-c$$$C Calculate the radial part of the gradient
-c$$$ gg(1)=xj*fac
-c$$$ gg(2)=yj*fac
-c$$$ gg(3)=zj*fac
-c$$$C Calculate angular part of the gradient.
-c$$$ call sc_grad
-c$$$c$$$ enddo ! j
-c$$$c$$$ enddo ! iint
-c$$$c$$$ enddo ! i
-c$$$ energy_dec=.false.
-c$$$ return
-c$$$ end
-c$$$
-c$$$C-----------------------------------------------------------------------------
-c$$$
-c$$$ subroutine perturb_side_chain(i,angle)
-c$$$ implicit none
-c$$$
-c$$$c Includes
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.CHAIN'
-c$$$ include 'COMMON.GEO'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.LOCAL'
-c$$$ include 'COMMON.IOUNITS'
-c$$$
-c$$$c External functions
-c$$$ external ran_number
-c$$$ double precision ran_number
-c$$$
-c$$$c Input arguments
-c$$$ integer i
-c$$$ double precision angle ! In degrees
-c$$$
-c$$$c Local variables
-c$$$ integer i_sc
-c$$$ double precision rad_ang,rand_v(3),length,cost,sint
-c$$$
-c$$$
-c$$$ i_sc=i+nres
-c$$$ rad_ang=angle*deg2rad
-c$$$
-c$$$ length=0.0
-c$$$ do while (length.lt.0.01)
-c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
-c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
-c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
-c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
-c$$$ + rand_v(3)*rand_v(3)
-c$$$ length=sqrt(length)
-c$$$ rand_v(1)=rand_v(1)/length
-c$$$ rand_v(2)=rand_v(2)/length
-c$$$ rand_v(3)=rand_v(3)/length
-c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
-c$$$ + rand_v(3)*dc_norm(3,i_sc)
-c$$$ length=1.0D0-cost*cost
-c$$$ if (length.lt.0.0D0) length=0.0D0
-c$$$ length=sqrt(length)
-c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
-c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
-c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
-c$$$ enddo
-c$$$ rand_v(1)=rand_v(1)/length
-c$$$ rand_v(2)=rand_v(2)/length
-c$$$ rand_v(3)=rand_v(3)/length
-c$$$
-c$$$ cost=dcos(rad_ang)
-c$$$ sint=dsin(rad_ang)
-c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
-c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
-c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
-c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
-c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
-c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
-c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
-c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
-c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
-c$$$
-c$$$ call chainbuild_cart
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-c$$$c----------------------------------------------------------------------------
-c$$$
-c$$$ subroutine ss_relax3(i_in,j_in)
-c$$$ implicit none
-c$$$
-c$$$c Includes
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.CHAIN'
-c$$$ include 'COMMON.IOUNITS'
-c$$$ include 'COMMON.INTERACT'
-c$$$
-c$$$c External functions
-c$$$ external ran_number
-c$$$ double precision ran_number
-c$$$
-c$$$c Input arguments
-c$$$ integer i_in,j_in
-c$$$
-c$$$c Local variables
-c$$$ double precision energy_sc(0:n_ene),etot
-c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
-c$$$ double precision ang_pert,rand_fact,exp_fact,beta
-c$$$ integer n,i_pert,i
-c$$$ logical notdone
-c$$$
-c$$$
-c$$$ beta=1.0D0
-c$$$
-c$$$ mask_r=.true.
-c$$$ do i=nnt,nct
-c$$$ mask_side(i)=0
-c$$$ enddo
-c$$$ mask_side(i_in)=1
-c$$$ mask_side(j_in)=1
-c$$$
-c$$$ call etotal_sc(energy_sc)
-c$$$ etot=energy_sc(0)
-c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
-c$$$c + energy_sc(1),energy_sc(12)
-c$$$
-c$$$ notdone=.true.
-c$$$ n=0
-c$$$ do while (notdone)
-c$$$ if (mod(n,2).eq.0) then
-c$$$ i_pert=i_in
-c$$$ else
-c$$$ i_pert=j_in
-c$$$ endif
-c$$$ n=n+1
-c$$$
-c$$$ do i=1,3
-c$$$ org_dc(i)=dc(i,i_pert+nres)
-c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
-c$$$ org_c(i)=c(i,i_pert+nres)
-c$$$ enddo
-c$$$ ang_pert=ran_number(0.0D0,3.0D0)
-c$$$ call perturb_side_chain(i_pert,ang_pert)
-c$$$ call etotal_sc(energy_sc)
-c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
-c$$$ rand_fact=ran_number(0.0D0,1.0D0)
-c$$$ if (rand_fact.lt.exp_fact) then
-c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
-c$$$c + energy_sc(1),energy_sc(12)
-c$$$ etot=energy_sc(0)
-c$$$ else
-c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
-c$$$c + energy_sc(1),energy_sc(12)
-c$$$ do i=1,3
-c$$$ dc(i,i_pert+nres)=org_dc(i)
-c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
-c$$$ c(i,i_pert+nres)=org_c(i)
-c$$$ enddo
-c$$$ endif
-c$$$
-c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
-c$$$ enddo
-c$$$
-c$$$ mask_r=.false.
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-c$$$c----------------------------------------------------------------------------
-c$$$
-c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
-c$$$ implicit none
-c$$$ include 'DIMENSIONS'
-c$$$ integer liv,lv
-c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
-c$$$*********************************************************************
-c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
-c$$$* the calling subprogram. *
-c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
-c$$$* calculated in the usual pythagorean way. *
-c$$$* absolute convergence occurs when the function is within v(31) of *
-c$$$* zero. unless you know the minimum value in advance, abs convg *
-c$$$* is probably not useful. *
-c$$$* relative convergence is when the model predicts that the function *
-c$$$* will decrease by less than v(32)*abs(fun). *
-c$$$*********************************************************************
-c$$$ include 'COMMON.IOUNITS'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.GEO'
-c$$$ include 'COMMON.MINIM'
-c$$$ include 'COMMON.CHAIN'
-c$$$
-c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
-c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
-c$$$ + orig_ss_dist(maxres2,maxres2)
-c$$$
-c$$$ double precision etot
-c$$$ integer iretcode,nfun,i_in,j_in
-c$$$
-c$$$ external dist
-c$$$ double precision dist
-c$$$ external ss_func,fdum
-c$$$ double precision ss_func,fdum
-c$$$
-c$$$ integer iv(liv),uiparm(2)
-c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
-c$$$ integer i,j,k
-c$$$
-c$$$
-c$$$ call deflt(2,iv,liv,lv,v)
-c$$$* 12 means fresh start, dont call deflt
-c$$$ iv(1)=12
-c$$$* max num of fun calls
-c$$$ if (maxfun.eq.0) maxfun=500
-c$$$ iv(17)=maxfun
-c$$$* max num of iterations
-c$$$ if (maxmin.eq.0) maxmin=1000
-c$$$ iv(18)=maxmin
-c$$$* controls output
-c$$$ iv(19)=2
-c$$$* selects output unit
-c$$$c iv(21)=iout
-c$$$ iv(21)=0
-c$$$* 1 means to print out result
-c$$$ iv(22)=0
-c$$$* 1 means to print out summary stats
-c$$$ iv(23)=0
-c$$$* 1 means to print initial x and d
-c$$$ iv(24)=0
-c$$$* min val for v(radfac) default is 0.1
-c$$$ v(24)=0.1D0
-c$$$* max val for v(radfac) default is 4.0
-c$$$ v(25)=2.0D0
-c$$$c v(25)=4.0D0
-c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
-c$$$* the sumsl default is 0.1
-c$$$ v(26)=0.1D0
-c$$$* false conv if (act fnctn decrease) .lt. v(34)
-c$$$* the sumsl default is 100*machep
-c$$$ v(34)=v(34)/100.0D0
-c$$$* absolute convergence
-c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
-c$$$ v(31)=tolf
-c$$$ v(31)=1.0D-1
-c$$$* relative convergence
-c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
-c$$$ v(32)=rtolf
-c$$$ v(32)=1.0D-1
-c$$$* controls initial step size
-c$$$ v(35)=1.0D-1
-c$$$* large vals of d correspond to small components of step
-c$$$ do i=1,6*nres
-c$$$ d(i)=1.0D0
-c$$$ enddo
-c$$$
-c$$$ do i=0,2*nres
-c$$$ do j=1,3
-c$$$ orig_ss_dc(j,i)=dc(j,i)
-c$$$ enddo
-c$$$ enddo
-c$$$ call geom_to_var(nvar,orig_ss_var)
-c$$$
-c$$$ do i=1,nres
-c$$$ do j=i,nres
-c$$$ orig_ss_dist(j,i)=dist(j,i)
-c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
-c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
-c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
-c$$$ enddo
-c$$$ enddo
-c$$$
-c$$$ k=0
-c$$$ do i=1,nres-1
-c$$$ do j=1,3
-c$$$ k=k+1
-c$$$ x(k)=dc(j,i)
-c$$$ enddo
-c$$$ enddo
-c$$$ do i=2,nres-1
-c$$$ if (ialph(i,1).gt.0) then
-c$$$ do j=1,3
-c$$$ k=k+1
-c$$$ x(k)=dc(j,i+nres)
-c$$$ enddo
-c$$$ endif
-c$$$ enddo
-c$$$
-c$$$ uiparm(1)=i_in
-c$$$ uiparm(2)=j_in
-c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
-c$$$ etot=v(10)
-c$$$ iretcode=iv(1)
-c$$$ nfun=iv(6)+iv(30)
-c$$$
-c$$$ k=0
-c$$$ do i=1,nres-1
-c$$$ do j=1,3
-c$$$ k=k+1
-c$$$ dc(j,i)=x(k)
-c$$$ enddo
-c$$$ enddo
-c$$$ do i=2,nres-1
-c$$$ if (ialph(i,1).gt.0) then
-c$$$ do j=1,3
-c$$$ k=k+1
-c$$$ dc(j,i+nres)=x(k)
-c$$$ enddo
-c$$$ endif
-c$$$ enddo
-c$$$ call chainbuild_cart
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-c$$$C-----------------------------------------------------------------------------
-c$$$
-c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
-c$$$ implicit none
-c$$$ include 'DIMENSIONS'
-c$$$ include 'COMMON.DERIV'
-c$$$ include 'COMMON.IOUNITS'
-c$$$ include 'COMMON.VAR'
-c$$$ include 'COMMON.CHAIN'
-c$$$ include 'COMMON.INTERACT'
-c$$$ include 'COMMON.SBRIDGE'
-c$$$
-c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
-c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
-c$$$ + orig_ss_dist(maxres2,maxres2)
-c$$$
-c$$$ integer n
-c$$$ double precision x(maxres6)
-c$$$ integer nf
-c$$$ double precision f
-c$$$ integer uiparm(2)
-c$$$ real*8 urparm(1)
-c$$$ external ufparm
-c$$$ double precision ufparm
-c$$$
-c$$$ external dist
-c$$$ double precision dist
-c$$$
-c$$$ integer i,j,k,ss_i,ss_j
-c$$$ double precision tempf,var(maxvar)
-c$$$
-c$$$
-c$$$ ss_i=uiparm(1)
-c$$$ ss_j=uiparm(2)
-c$$$ f=0.0D0
-c$$$
-c$$$ k=0
-c$$$ do i=1,nres-1
-c$$$ do j=1,3
-c$$$ k=k+1
-c$$$ dc(j,i)=x(k)
-c$$$ enddo
-c$$$ enddo
-c$$$ do i=2,nres-1
-c$$$ if (ialph(i,1).gt.0) then
-c$$$ do j=1,3
-c$$$ k=k+1
-c$$$ dc(j,i+nres)=x(k)
-c$$$ enddo
-c$$$ endif
-c$$$ enddo
-c$$$ call chainbuild_cart
-c$$$
-c$$$ call geom_to_var(nvar,var)
-c$$$
-c$$$c Constraints on all angles
-c$$$ do i=1,nvar
-c$$$ tempf=var(i)-orig_ss_var(i)
-c$$$ f=f+tempf*tempf
-c$$$ enddo
-c$$$
-c$$$c Constraints on all distances
-c$$$ do i=1,nres-1
-c$$$ if (i.gt.1) then
-c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
-c$$$ f=f+tempf*tempf
-c$$$ endif
-c$$$ do j=i+1,nres
-c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
-c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
-c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
-c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
-c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
-c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
-c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
-c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
-c$$$ enddo
-c$$$ enddo
-c$$$
-c$$$c Constraints for the relevant CYS-CYS
-c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
-c$$$ f=f+tempf*tempf
-c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
-c$$$
-c$$$c$$$ if (nf.ne.nfl) then
-c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
-c$$$c$$$ + f,dist(5+nres,14+nres)
-c$$$c$$$ endif
-c$$$
-c$$$ nfl=nf
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-c$$$C-----------------------------------------------------------------------------
-c$$$C-----------------------------------------------------------------------------
- subroutine triple_ssbond_ene(resi,resj,resk,eij)