rescode.f
setup_var.f
srtclust.f
+ ssMD.F
timing.F
track.F
wrtclust.f
probabl.F
read_coords.F
readrtns.F
+ ssMD.F
timing.F
track.F
work_partition.F
#=========================================
# Additional flags
#=========================================
-set(CPPFLAGS "${CPPFLAGS} -DUNRES -DISNAN" )
+set(CPPFLAGS "${CPPFLAGS} -DUNRES -DISNAN -DCLUST" )
#=========================================
# Compiler specific flags
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*fact(1)*ees+wvdwpp*evdw1
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
#else
etot=wsc*evdw+wscp*evdw2+welec*fact(1)*(ees+evdw1)
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
energia(18)=estr
energia(19)=esccor
energia(20)=edihcnstr
+cc if (dyn_ss) call dyn_set_nss
c detecting NaNQ
i=0
#ifdef WINPGI
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SBRIDGE'
logical lprn
common /srutu/icall
integer icant
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ call dyn_ssbond_ene(i,j,evdwij)
+ evdw=evdw+evdwij
+c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+c & 'evdw',i,j,evdwij,' ss'
+ ELSE
ind=ind+1
itypj=itype(j)
dscj_inv=vbld_inv(j+nres)
C Calculate angular part of the gradient.
call sc_grad
endif
+ ENDIF ! SSBOND
enddo ! j
enddo ! iint
enddo ! i
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SBRIDGE'
common /srutu/ icall
logical lprn
integer icant
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
+C in case of diagnostics write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j)
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ call dyn_ssbond_ene(i,j,evdwij)
+ evdw=evdw+evdwij
+c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+c & 'evdw',i,j,evdwij,' ss'
+ ELSE
ind=ind+1
itypj=itype(j)
dscj_inv=vbld_inv(j+nres)
C Calculate angular part of the gradient.
call sc_grad
endif
+ ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
+ if (.not.dyn_ss .and. i.le.nss) then
+C 15/02/13 CC dynamic SSbond - additional check
if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
cd write (iout,*) "eij",eij
+ endif
else if (ii.gt.nres .and. jj.gt.nres) then
c Restraints from contact prediction
dd=dist(ii,jj)
deltat12=om2-om1+2.0d0
cosphi=om12-om1*om2
eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
- & +akct*deltad*deltat12
+ & +akct*deltad*deltat12+ebr
& +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
C Set lprn=.true. for debugging
lprn=.false.
eturn6=0.0d0
+ ecorr6=0.0d0
#ifdef MPL
n_corr=0
n_corr1=0
cd write(2,*)'ijkl',i,j,i+1,j1
if (wcorr6.gt.0.0d0 .and. (j.ne.i+4 .or. j1.ne.i+3
& .or. wturn6.eq.0.0d0))then
-cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
- ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
-cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
-cd & 'ecorr6=',ecorr6
+c write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
+c ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
+c write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
+c & 'ecorr6=',ecorr6, wcorr6
cd write (iout,'(4e15.5)') sred_geom,
cd & dabs(eello4(i,j,i+1,j1,jj,kk)),
cd & dabs(eello5(i,j,i+1,j1,jj,kk)),
- double precision ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,dhpb,
- & dhpb1,forcon,weidis
- integer ns,nss,nfree,iss,ihpb,jhpb,nhpb,link_start,link_end,
- & ibecarb
- common /sbridge/ ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,ns,nss,
- & nfree,iss(maxss)
+ double precision ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss
+ integer ns,nss,nfree,iss
+ common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,
+ & ns,nss,nfree,iss(maxss)
+ double precision dhpb,dhpb1,forcon
+ integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb
common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
& ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb
+ double precision weidis
common /restraints/ weidis
+ integer link_start,link_end
common /links_split/ link_start,link_end
+ double precision Ht,dyn_ssbond_ij
+ logical dyn_ss,dyn_ss_mask
+ common /dyn_ssbond/ dyn_ssbond_ij(maxres,maxres),
+ & idssb(maxdim),jdssb(maxdim),
+ & Ht,dyn_ss,dyn_ss_mask(maxres)
ihpb(i)=0
jhpb(i)=0
enddo
+ do i=1,maxres
+ dyn_ss_mask(i)=.false.
+ enddo
C
C Initialize timing.
C
do ii=1,nss
if (ihpb(ii).eq.i+nres) then
scheck=.true.
+ if (dyn_ss) go to 10
jj=jhpb(ii)-nres
goto 10
endif
#ifdef MPI
call work_partition(.true.,ncon)
#endif
-
call probabl(iT,ncon_work,ncon,*20)
if (ncon_work.lt.2) then
C
C Define the constants of the disulfide bridge
C
- ebr=-5.50D0
+C ebr=-5.50D0
c
c Old arbitrary potential - commented out.
c
c energy surface of diethyl disulfide.
c A. Liwo and U. Kozlowska, 11/24/03
c
- D0CM = 3.78d0
- AKCM = 15.1d0
- AKTH = 11.0d0
- AKCT = 12.0d0
- V1SS =-1.08d0
- V2SS = 7.61d0
- V3SS = 13.7d0
-
- write (iout,'(/a)') "Disulfide bridge parameters:"
- write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
- write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
- write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
- write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
- & ' v3ss:',v3ss
+C D0CM = 3.78d0
+C AKCM = 15.1d0
+C AKTH = 11.0d0
+C AKCT = 12.0d0
+C V1SS =-1.08d0
+C V2SS = 7.61d0
+C V3SS = 13.7d0
return
end
c do i=1,ncon
c write (iout,*) i,list_conf(i)
c enddo
+c do i=1,ncon
+c write(iout,*) "entrop before", entfac(i),i
+c enddo
+
#ifdef MPI
write (iout,*) me," indstart",indstart(me)," indend",indend(me)
call daread_ccoords(indstart(me),indend(me))
#endif
+c do i=1,ncon
+c write(iout,*) "entrop after", entfac(i),i
+c enddo
+
c write (iout,*) "ncon",ncon
+
temper=1.0d0/(beta_h(ib)*1.987D-3)
c write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
c quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
call int_from_cart1(.false.)
call etotal(energia(0),fT)
totfree(i)=energia(0)
+c#define DEBUG
#ifdef DEBUG
- write (iout,*) i," energia",(energia(j),j=0,21)
+ write (iout,*) i," energia",(energia(j),j=0,20)
call enerprint(energia(0),ft)
call flush(iout)
#endif
+c#undef DEBUG
do k=1,max_ene
enetb(k,i)=energia(k)
enddo
ecorr=enetb(4,i)
ecorr5=enetb(5,i)
ecorr6=enetb(6,i)
+cc if (wcorr6.eq.0) ecorr6=0.0d0
eel_loc=enetb(7,i)
eello_turn3=enetb(8,i)
eello_turn4=enetb(9,i)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+ft(1)*welec*ees+wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
#else
etot=wsc*evdw+wscp*evdw2+ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
Fdimless(i)=beta_h(ib)*etot+entfac(ii)
totfree(i)=etot
#ifdef DEBUG
- write (iout,*) i,ii,ib,
+
+ write (iout,*) "etrop", i,ii,ib,
& 1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
& entfac(ii),Fdimless(i)
#endif
& MPI_COMM_WORLD, IERROR)
if (me.eq.Master) then
#endif
+c#define DEBUG
#ifdef DEBUG
write (iout,*) "The FDIMLESS array before sorting"
do i=1,ncon
write (iout,*) i,list_conf(i),fdimless(i)
enddo
#endif
+c#undef DEBUG
do i=1,ncon
totfree(i)=fdimless(i)
enddo
qfree=0.0d0
do i=1,ncon
- qfree=qfree+exp(-fdimless(i)+fdimless(1))
+ qfree=qfree+dexp(dble(-fdimless(i)+fdimless(1)))
enddo
c write (iout,*) "qfree",qfree
nlist=1
sumprob=0.0
do i=1,min0(ncon,maxstr_proc)-1
- sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
+ sumprob=sumprob+dexp(dble(-fdimless(i)+fdimless(1)))/qfree
+c#define DEBUG
#ifdef DEBUG
- write (iout,*) i,ib,beta_h(ib),
+ write (iout,*) 'i=',i,ib,beta_h(ib),
& 1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
& totfree(list_conf(i)),
& -entfac(list_conf(i)),fdimless(i),sumprob
#endif
+c#undef DEBUG
if (sumprob.gt.prob_limit) goto 122
c if (sumprob.gt.1.00d0) goto 122
nlist=nlist+1
call xdrfint_(ixdrf, nss, iret)
if (iret.eq.0) goto 101
do j=1,nss
- call xdrfint_(ixdrf, ihpb(j), iret)
- if (iret.eq.0) goto 101
- call xdrfint_(ixdrf, jhpb(j), iret)
- if (iret.eq.0) goto 101
+cc if (dyn_ss) then
+cc call xdrfint_(ixdrf, idssb(j), iret)
+cc if (iret.eq.0) goto 101
+cc call xdrfint_(ixdrf, jdssb(j), iret)
+cc if (iret.eq.0) goto 101
+cc idssb(j)=idssb(j)-nres
+cc jdssb(j)=jdssb(j)-nres
+cc else
+ call xdrfint_(ixdrf, ihpb(j), iret)
+ if (iret.eq.0) goto 101
+ call xdrfint_(ixdrf, jhpb(j), iret)
+ if (iret.eq.0) goto 101
+cc endif
enddo
call xdrffloat_(ixdrf,reini,iret)
if (iret.eq.0) goto 101
call flush(iout)
if (iret.eq.0) goto 101
do k=1,nss
- call xdrfint(ixdrf, ihpb(k), iret)
- if (iret.eq.0) goto 101
- call xdrfint(ixdrf, jhpb(k), iret)
- if (iret.eq.0) goto 101
+cc if (dyn_ss) then
+cc call xdrfint(ixdrf, idssb(k), iret)
+cc if (iret.eq.0) goto 101
+cc call xdrfint(ixdrf, jdssb(k), iret)
+cc if (iret.eq.0) goto 101
+cc idssb(k)=idssb(k)-nres
+cc jdssb(k)=jdssb(k)-nres
+cc write(iout,*) "TUTU", idssb(k),jdssb(k)
+cc else
+ call xdrfint(ixdrf, ihpb(k), iret)
+ if (iret.eq.0) goto 101
+ call xdrfint(ixdrf, jhpb(k), iret)
+ if (iret.eq.0) goto 101
+cc endif
enddo
call xdrffloat(ixdrf,reini,iret)
if (iret.eq.0) goto 101
if (iret.eq.0) goto 101
#endif
energy(jj+1)=reini
- entfac(jj+1)=refree
+cc write(iout,*) 'reini=', reini, jj+1
+ entfac(jj+1)=dble(refree)
+cc write(iout,*) 'refree=', refree,jj+1
rmstb(jj+1)=rmsdev
do k=1,nres
do l=1,3
write (iout,*) "Reading binary file, record",iii," ii",ii
call flush(iout)
#endif
+ if (dyn_ss) then
+ read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),
+ & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),
+c & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),
+ & entfac(ii),rmstb(ii)
+ else
read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),
& ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),
& nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),
& entfac(ii),rmstb(ii)
+ endif
#ifdef DEBUG
write (iout,*) ii,iii,ij,entfac(ii)
write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
write (iout,*) "Writing binary file, record",iii," ii",ii
call flush(iout)
#endif
+ if (dyn_ss) then
+ write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),
+ & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),
+c & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij))
+ & entfac(ii),rmstb(ii)
+ else
write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),
& ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),
& nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),
& entfac(ii),rmstb(ii)
+ endif
#ifdef DEBUG
write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,
refstr=(index(controlcard,'REFSTR').gt.0)
write (iout,*) "REFSTR",refstr
pdbref=(index(controlcard,'PDBREF').gt.0)
+ dyn_ss=(index(controlcard,'DYN_SS').gt.0)
iscode=index(controlcard,'ONE_LETTER')
tree=(index(controlcard,'MAKE_TREE').gt.0)
min_var=(index(controlcard,'MINVAR').gt.0)
call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
call reada(weightcard,'WSCCOR',wsccor,1.0D0)
+ call reada(weightcard,"D0CM",d0cm,3.78d0)
+ call reada(weightcard,"AKCM",akcm,15.1d0)
+ call reada(weightcard,"AKTH",akth,11.0d0)
+ call reada(weightcard,"AKCT",akct,12.0d0)
+ call reada(weightcard,"V1SS",v1ss,-1.08d0)
+ call reada(weightcard,"V2SS",v2ss,7.61d0)
+ call reada(weightcard,"V3SS",v3ss,13.7d0)
+ call reada(weightcard,"EBR",ebr,-5.50D0)
if (index(weightcard,'SOFT').gt.0) ipot=6
C 12/1/95 Added weight for the multi-body term WCORR
call reada(weightcard,'WCORRH',wcorr,1.0D0)
+ do i=1,maxres-1
+ do j=i+1,maxres
+ dyn_ssbond_ij(i,j)=1.0d300
+ enddo
+ enddo
+ call reada(weightcard,"HT",Ht,0.0D0)
+ if (dyn_ss) then
+ ss_depth=ebr/wsc-0.25*eps(1,1)
+ Ht=Ht/wsc-0.25*eps(1,1)
+ akcm=akcm*wstrain/wsc
+ akth=akth*wstrain/wsc
+ akct=akct*wstrain/wsc
+ v1ss=v1ss*wstrain/wsc
+ v2ss=v2ss*wstrain/wsc
+ v3ss=v3ss*wstrain/wsc
+ else
+ ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+ endif
+ write (iout,'(/a)') "Disulfide bridge parameters:"
+ write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
+ write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
+ write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
+ write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
+ write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
+ & ' v3ss:',v3ss
+ write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
if (wcorr4.gt.0.0d0) wcorr=wcorr4
weights(1)=wsc
weights(2)=wscp
enddo
endif
endif
+ if (ns.gt.0.and.dyn_ss) then
+ do i=nss+1,nhpb
+ ihpb(i-nss)=ihpb(i)
+ jhpb(i-nss)=jhpb(i)
+ forcon(i-nss)=forcon(i)
+ dhpb(i-nss)=dhpb(i)
+ enddo
+ nhpb=nhpb-nss
+ nss=0
+ call hpb_partition
+ do i=1,ns
+ dyn_ss_mask(iss(i))=.true.
+c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
+ enddo
+ endif
+ print *, "Leaving brigde read"
return
end
c----------------------------------------------------------------------------
--- /dev/null
+c----------------------------------------------------------------------------
+ subroutine check_energies
+c implicit none
+
+c Includes
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.LOCAL'
+ include 'COMMON.GEO'
+
+c External functions
+ double precision ran_number
+ external ran_number
+
+c Local variables
+ integer i,j,k,l,lmax,p,pmax
+ double precision rmin,rmax
+ double precision eij
+
+ double precision d
+ double precision wi,rij,tj,pj
+
+
+c return
+
+ i=5
+ j=14
+
+ d=dsc(1)
+ rmin=2.0D0
+ rmax=12.0D0
+
+ lmax=10000
+ pmax=1
+
+ do k=1,3
+ c(k,i)=0.0D0
+ c(k,j)=0.0D0
+ c(k,nres+i)=0.0D0
+ c(k,nres+j)=0.0D0
+ enddo
+
+ do l=1,lmax
+
+ct wi=ran_number(0.0D0,pi)
+c wi=ran_number(0.0D0,pi/6.0D0)
+c wi=0.0D0
+ct tj=ran_number(0.0D0,pi)
+ct pj=ran_number(0.0D0,pi)
+c pj=ran_number(0.0D0,pi/6.0D0)
+c pj=0.0D0
+
+ do p=1,pmax
+ct rij=ran_number(rmin,rmax)
+
+ c(1,j)=d*sin(pj)*cos(tj)
+ c(2,j)=d*sin(pj)*sin(tj)
+ c(3,j)=d*cos(pj)
+
+ c(3,nres+i)=-rij
+
+ c(1,i)=d*sin(wi)
+ c(3,i)=-rij-d*cos(wi)
+
+ do k=1,3
+ dc(k,nres+i)=c(k,nres+i)-c(k,i)
+ dc_norm(k,nres+i)=dc(k,nres+i)/d
+ dc(k,nres+j)=c(k,nres+j)-c(k,j)
+ dc_norm(k,nres+j)=dc(k,nres+j)/d
+ enddo
+
+ call dyn_ssbond_ene(i,j,eij)
+ enddo
+ enddo
+
+ call exit(1)
+
+ return
+ end
+
+C-----------------------------------------------------------------------------
+
+ subroutine dyn_ssbond_ene(resi,resj,eij)
+c implicit none
+
+c Includes
+ include 'DIMENSIONS'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+ include 'COMMON.MD'
+#endif
+#endif
+
+c External functions
+ double precision h_base
+ external h_base
+
+c Input arguments
+ integer resi,resj
+
+c Output arguments
+ double precision eij
+
+c Local variables
+ logical havebond
+c integer itypi,itypj,k,l
+ double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
+ double precision sig0ij,ljd,sig,fac,e1,e2
+ double precision dcosom1(3),dcosom2(3),ed
+ double precision pom1,pom2
+ double precision ljA,ljB,ljXs
+ double precision d_ljB(1:3)
+ double precision ssA,ssB,ssC,ssXs
+ double precision ssxm,ljxm,ssm,ljm
+ double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+ double precision f1,f2,h1,h2,hd1,hd2
+ double precision omega,delta_inv,deltasq_inv,fac1,fac2
+c-------FIRST METHOD
+ double precision xm,d_xm(1:3)
+c-------END FIRST METHOD
+c-------SECOND METHOD
+c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
+c-------END SECOND METHOD
+
+c-------TESTING CODE
+ logical checkstop,transgrad
+ common /sschecks/ checkstop,transgrad
+
+ integer icheck,nicheck,jcheck,njcheck
+ double precision echeck(-1:1),deps,ssx0,ljx0
+c-------END TESTING CODE
+
+
+ i=resi
+ j=resj
+
+ itypi=itype(i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+
+ itypj=itype(j)
+ xj=c(1,nres+j)-c(1,nres+i)
+ yj=c(2,nres+j)-c(2,nres+i)
+ zj=c(3,nres+j)-c(3,nres+i)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ dscj_inv=vbld_inv(j+nres)
+
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
+c The following are set in sc_angular
+c erij(1)=xj*rij
+c erij(2)=yj*rij
+c erij(3)=zj*rij
+c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+c om12=dxi*dxj+dyi*dyj+dzi*dzj
+ call sc_angular
+ rij=1.0D0/rij ! Reset this so it makes sense
+
+ sig0ij=sigma(itypi,itypj)
+ sig=sig0ij*dsqrt(1.0D0/sigsq)
+
+ ljXs=sig-sig0ij
+ ljA=eps1*eps2rt**2*eps3rt**2
+ ljB=ljA*bb(itypi,itypj)
+ ljA=ljA*aa(itypi,itypj)
+ ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+
+ ssXs=d0cm
+ deltat1=1.0d0-om1
+ deltat2=1.0d0+om2
+ deltat12=om2-om1+2.0d0
+ cosphi=om12-om1*om2
+ ssA=akcm
+ ssB=akct*deltat12
+ ssC=ss_depth
+ & +akth*(deltat1*deltat1+deltat2*deltat2)
+ & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
+ ssxm=ssXs-0.5D0*ssB/ssA
+
+c-------TESTING CODE
+c$$$c Some extra output
+c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
+c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
+c$$$ if (ssx0.gt.0.0d0) then
+c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
+c$$$ else
+c$$$ ssx0=ssxm
+c$$$ endif
+c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
+c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
+c$$$ return
+c-------END TESTING CODE
+
+c-------TESTING CODE
+c Stop and plot energy and derivative as a function of distance
+ if (checkstop) then
+ ssm=ssC-0.25D0*ssB*ssB/ssA
+ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+ if (ssm.lt.ljm .and.
+ & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
+ nicheck=1000
+ njcheck=1
+ deps=0.5d-7
+ else
+ checkstop=.false.
+ endif
+ endif
+ if (.not.checkstop) then
+ nicheck=0
+ njcheck=-1
+ endif
+
+ do icheck=0,nicheck
+ do jcheck=-1,njcheck
+ if (checkstop) rij=(ssxm-1.0d0)+
+ & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
+c-------END TESTING CODE
+
+ if (rij.gt.ljxm) then
+ havebond=.false.
+ ljd=rij-ljXs
+ fac=(1.0D0/ljd)**expon
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ eij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=eij*eps3rt
+ eps3der=eij*eps2rt
+ eij=eij*eps2rt*eps3rt
+
+ sigder=-sig/sigsq
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ ed=-expon*(e1+eij)/ljd
+ sigder=ed*sigder
+ eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
+ eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
+ eom12=eij*eps1_om12+eps2der*eps2rt_om12
+ & -2.0D0*alf12*eps3der+sigder*sigsq_om12
+ else if (rij.lt.ssxm) then
+ havebond=.true.
+ ssd=rij-ssXs
+ eij=ssA*ssd*ssd+ssB*ssd+ssC
+
+ ed=2*akcm*ssd+akct*deltat12
+ pom1=akct*ssd
+ pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
+ eom1=-2*akth*deltat1-pom1-om2*pom2
+ eom2= 2*akth*deltat2+pom1-om1*pom2
+ eom12=pom2
+ else
+ omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
+
+ d_ssxm(1)=0.5D0*akct/ssA
+ d_ssxm(2)=-d_ssxm(1)
+ d_ssxm(3)=0.0D0
+
+ d_ljxm(1)=sig0ij/sqrt(sigsq**3)
+ d_ljxm(2)=d_ljxm(1)*sigsq_om2
+ d_ljxm(3)=d_ljxm(1)*sigsq_om12
+ d_ljxm(1)=d_ljxm(1)*sigsq_om1
+
+c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
+ xm=0.5d0*(ssxm+ljxm)
+ do k=1,3
+ d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
+ enddo
+ if (rij.lt.xm) then
+ havebond=.true.
+ ssm=ssC-0.25D0*ssB*ssB/ssA
+ d_ssm(1)=0.5D0*akct*ssB/ssA
+ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
+ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
+ d_ssm(3)=omega
+ f1=(rij-xm)/(ssxm-xm)
+ f2=(rij-ssxm)/(xm-ssxm)
+ h1=h_base(f1,hd1)
+ h2=h_base(f2,hd2)
+ eij=ssm*h1+Ht*h2
+ delta_inv=1.0d0/(xm-ssxm)
+ deltasq_inv=delta_inv*delta_inv
+ fac=ssm*hd1-Ht*hd2
+ fac1=deltasq_inv*fac*(xm-rij)
+ fac2=deltasq_inv*fac*(rij-ssxm)
+ ed=delta_inv*(Ht*hd2-ssm*hd1)
+ eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
+ eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
+ eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
+ else
+ havebond=.false.
+ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+ d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
+ d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
+ d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
+ + alf12/eps3rt)
+ d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
+ f1=(rij-ljxm)/(xm-ljxm)
+ f2=(rij-xm)/(ljxm-xm)
+ h1=h_base(f1,hd1)
+ h2=h_base(f2,hd2)
+ eij=Ht*h1+ljm*h2
+ delta_inv=1.0d0/(ljxm-xm)
+ deltasq_inv=delta_inv*delta_inv
+ fac=Ht*hd1-ljm*hd2
+ fac1=deltasq_inv*fac*(ljxm-rij)
+ fac2=deltasq_inv*fac*(rij-xm)
+ ed=delta_inv*(ljm*hd2-Ht*hd1)
+ eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
+ eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
+ eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
+ endif
+c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
+
+c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
+c$$$ ssd=rij-ssXs
+c$$$ ljd=rij-ljXs
+c$$$ fac1=rij-ljxm
+c$$$ fac2=rij-ssxm
+c$$$
+c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
+c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
+c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
+c$$$
+c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
+c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
+c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
+c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
+c$$$ d_ssm(3)=omega
+c$$$
+c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
+c$$$ do k=1,3
+c$$$ d_ljm(k)=ljm*d_ljB(k)
+c$$$ enddo
+c$$$ ljm=ljm*ljB
+c$$$
+c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
+c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
+c$$$ d_ss(2)=akct*ssd
+c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
+c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
+c$$$ d_ss(3)=omega
+c$$$
+c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
+c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
+c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
+c$$$ do k=1,3
+c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
+c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
+c$$$ enddo
+c$$$ ljf=ljm+ljf*ljB*fac1*fac1
+c$$$
+c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
+c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
+c$$$ h1=h_base(f1,hd1)
+c$$$ h2=h_base(f2,hd2)
+c$$$ eij=ss*h1+ljf*h2
+c$$$ delta_inv=1.0d0/(ljxm-ssxm)
+c$$$ deltasq_inv=delta_inv*delta_inv
+c$$$ fac=ljf*hd2-ss*hd1
+c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
+c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
+c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
+c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
+c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
+c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
+c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
+c$$$
+c$$$ havebond=.false.
+c$$$ if (ed.gt.0.0d0) havebond=.true.
+c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
+
+ endif
+
+ if (havebond) then
+#ifndef CLUST
+#ifndef WHAM
+c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
+c write(iout,'(a15,f12.2,f8.1,2i5)')
+c & "SSBOND_E_FORM",totT,t_bath,i,j
+c endif
+#endif
+#endif
+ write(iout,*), 'DYN_SS_BOND',i,j,eij
+ dyn_ssbond_ij(i,j)=eij
+ else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
+ dyn_ssbond_ij(i,j)=1.0d300
+#ifndef CLUST
+#ifndef WHAM
+c write(iout,'(a15,f12.2,f8.1,2i5)')
+c & "SSBOND_E_BREAK",totT,t_bath,i,j
+#endif
+#endif
+ endif
+
+c-------TESTING CODE
+ if (checkstop) then
+ if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
+ & "CHECKSTOP",rij,eij,ed
+ echeck(jcheck)=eij
+ endif
+ enddo
+ if (checkstop) then
+ write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
+ endif
+ enddo
+ if (checkstop) then
+ transgrad=.true.
+ checkstop=.false.
+ endif
+c-------END TESTING CODE
+
+ do k=1,3
+ dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
+ dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
+ enddo
+ do k=1,3
+ gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ enddo
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
+ & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
+ & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+cgrad do k=i,j-1
+cgrad do l=1,3
+cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
+cgrad enddo
+cgrad enddo
+
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+
+ return
+ end
+
+C-----------------------------------------------------------------------------
+
+ double precision function h_base(x,deriv)
+c A smooth function going 0->1 in range [0,1]
+c It should NOT be called outside range [0,1], it will not work there.
+ implicit none
+
+c Input arguments
+ double precision x
+
+c Output arguments
+ double precision deriv
+
+c Local variables
+ double precision xsq
+
+
+c Two parabolas put together. First derivative zero at extrema
+c$$$ if (x.lt.0.5D0) then
+c$$$ h_base=2.0D0*x*x
+c$$$ deriv=4.0D0*x
+c$$$ else
+c$$$ deriv=1.0D0-x
+c$$$ h_base=1.0D0-2.0D0*deriv*deriv
+c$$$ deriv=4.0D0*deriv
+c$$$ endif
+
+c Third degree polynomial. First derivative zero at extrema
+ h_base=x*x*(3.0d0-2.0d0*x)
+ deriv=6.0d0*x*(1.0d0-x)
+
+c Fifth degree polynomial. First and second derivatives zero at extrema
+c$$$ xsq=x*x
+c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
+c$$$ deriv=x-1.0d0
+c$$$ deriv=deriv*deriv
+c$$$ deriv=30.0d0*xsq*deriv
+
+ return
+ end
+
+c----------------------------------------------------------------------------
+
+ subroutine dyn_set_nss
+c Adjust nss and other relevant variables based on dyn_ssbond_ij
+c implicit none
+
+c Includes
+ include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+#endif
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+#ifndef CLUST
+ integer maxprocs
+ parameter (maxprocs=128)
+#endif
+c#ifdef WHAM
+ integer max_fg_procs
+ parameter (max_fg_procs=128)
+c#endif
+cc include 'COMMON.SETUP'
+#ifndef CLUST
+#ifndef WHAM
+ include 'COMMON.MD'
+ include 'COMMON.SETUP'
+#endif
+#endif
+
+c Local variables
+ double precision emin
+ integer i,j,imin
+ integer diff,allflag(maxdim),allnss,
+ & allihpb(maxdim),alljhpb(maxdim),
+ & newnss,newihpb(maxdim),newjhpb(maxdim)
+ logical found
+ integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
+ integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
+
+ allnss=0
+ do i=1,nres-1
+ do j=i+1,nres
+ if (dyn_ssbond_ij(i,j).lt.1.0d300) then
+ allnss=allnss+1
+ allflag(allnss)=0
+ allihpb(allnss)=i
+ alljhpb(allnss)=j
+ endif
+ enddo
+ enddo
+
+cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
+
+ 1 emin=1.0d300
+ do i=1,allnss
+ if (allflag(i).eq.0 .and.
+ & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
+ emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
+ imin=i
+ endif
+ enddo
+ if (emin.lt.1.0d300) then
+ allflag(imin)=1
+ do i=1,allnss
+ if (allflag(i).eq.0 .and.
+ & (allihpb(i).eq.allihpb(imin) .or.
+ & alljhpb(i).eq.allihpb(imin) .or.
+ & allihpb(i).eq.alljhpb(imin) .or.
+ & alljhpb(i).eq.alljhpb(imin))) then
+ allflag(i)=-1
+ endif
+ enddo
+ goto 1
+ endif
+
+cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
+
+ newnss=0
+ do i=1,allnss
+ if (allflag(i).eq.1) then
+ newnss=newnss+1
+ newihpb(newnss)=allihpb(i)
+ newjhpb(newnss)=alljhpb(i)
+ endif
+ enddo
+
+#ifdef MPI
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(newnss,g_newnss,1,
+ & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+ call MPI_Gather(newnss,1,MPI_INTEGER,
+ & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_newnss(i-1)+displ(i-1)
+ enddo
+ call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
+ & g_newihpb,i_newnss,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
+ & g_newjhpb,i_newnss,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ if(fg_rank.eq.0) then
+c print *,'g_newnss',g_newnss
+c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
+c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
+ newnss=g_newnss
+ do i=1,newnss
+ newihpb(i)=g_newihpb(i)
+ newjhpb(i)=g_newjhpb(i)
+ enddo
+ endif
+ endif
+#endif
+
+ diff=newnss-nss
+
+cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
+
+ do i=1,nss
+ found=.false.
+ do j=1,newnss
+ if (idssb(i).eq.newihpb(j) .and.
+ & jdssb(i).eq.newjhpb(j)) found=.true.
+ enddo
+#ifndef CLUST
+#ifndef WHAM
+ if (.not.found.and.fg_rank.eq.0)
+ & write(iout,'(a15,f12.2,f8.1,2i5)')
+ & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
+#endif
+#endif
+ enddo
+
+ do i=1,newnss
+ found=.false.
+ do j=1,nss
+ if (newihpb(i).eq.idssb(j) .and.
+ & newjhpb(i).eq.jdssb(j)) found=.true.
+ enddo
+#ifndef CLUST
+#ifndef WHAM
+ if (.not.found.and.fg_rank.eq.0)
+ & write(iout,'(a15,f12.2,f8.1,2i5)')
+ & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
+#endif
+#endif
+ enddo
+
+ nss=newnss
+ do i=1,nss
+ idssb(i)=newihpb(i)
+ jdssb(i)=newjhpb(i)
+ enddo
+
+ return
+ end
+
+c----------------------------------------------------------------------------
+#ifdef NIEWIADOM
+c#ifdef WHAM
+ subroutine read_ssHist
+ implicit none
+
+c Includes
+ include 'DIMENSIONS'
+ include "DIMENSIONS.FREE"
+ include 'COMMON.FREE'
+
+c Local variables
+ integer i,j
+ character*80 controlcard
+
+ do i=1,dyn_nssHist
+ call card_concat(controlcard,.true.)
+ read(controlcard,*)
+ & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
+ enddo
+
+ return
+ end
+#endif
+
+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-----------------------------------------------------------------------------
rescode.f
setup_var.f
slices.F
+ ssMD.F
store_parm.F
timing.F
wham_calc1.F
readrtns_compar.F
readrtns.F
slices.F
+ ssMD.F
store_parm.F
timing.F
wham_calc1.F
#=========================================
# Additional flags
#=========================================
-set(CPPFLAGS "${CPPFLAGS} -DUNRES -DISNAN")
+set(CPPFLAGS "${CPPFLAGS} -DUNRES -DISNAN -DWHAM ")
#=========================================
& sigma_all(ntyp,ntyp,max_parm),r0_all(ntyp,ntyp,max_parm),
& chi_all(ntyp,ntyp,max_parm),chip_all(ntyp,max_parm),
& alp_all(ntyp,max_parm),ebr_all(max_parm),d0cm_all(max_parm),
+ & ss_depth_all(max_parm),ht_all(max_parm),
& akcm_all(max_parm),akth_all(max_parm),akct_all(max_parm),
& v1ss_all(max_parm),v2ss_all(max_parm),v3ss_all(max_parm),
& v1sccor_all(maxterm_sccor,3,ntyp,ntyp,max_parm),
& ctilde_all,dtilde_all,b1tilde_all,app_all,bpp_all,ael6_all,
& ael3_all,aad_all,bad_all,aa_all,bb_all,augm_all,
& eps_all,sigma_all,r0_all,chi_all,chip_all,alp_all,ebr_all,
+ & ss_depth_all,ht_all,
& d0cm_all,akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all,
& v1sccor_all,v2sccor_all,nbondterm_all,
& nlob_all,nlor_all,nterm_all,ntermd1_all,ntermd2_all,
********************************************************************************
c implicit real*8 (a-h,o-z)
C Max. number of processors.
-c parameter (maxprocs=128)
+C parameter (maxprocs=128)
C Max. number of fine-grain processors
-c parameter (max_fg_procs=maxprocs)
+C parameter (max_fg_procs=maxprocs)
C Max. number of coarse-grain processors
c parameter (max_cg_procs=maxprocs)
C Max. number of AA residues
real*4 csingle(3,maxres2)
character*64 nazwa,bprotfile_temp
character*3 liczba
- integer i,is,ie,j,ii,jj,k,kk,l,ll,mm,if
+ integer i,is,ie,j,ii,jj,k,kk,l,ll,mm,if,m
integer nrec,nlines,iscor,islice
double precision energ
integer ilen,iroof
& eini,efree,rmsdev,(prop(j),j=1,nQ),iscor
ii=ii+1
kk=kk+1
+ write(iout,*) 'BXWEJ',eini,l
+ flush(iout)
if (mod(kk,isampl(iparm)).eq.0) then
jj=jj+1
write(ientout,rec=jj)
& nss,(ihpb(k),jhpb(k),k=1,nss),
& eini,efree,rmsdev,(prop(j),j=1,nQ),iR,ib,iparm
#ifdef DEBUG
- do i=1,2*nres
+ do l=1,2*nres
do j=1,3
- c(j,i)=csingle(j,i)
+ c(j,l)=csingle(j,l)
enddo
enddo
call int_from_cart1(.false.)
write (iout,*) "Writing conformation, record",jj
write (iout,*) "Cartesian coordinates"
- write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
- write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
+ write (iout,'(8f10.5)') ((c(j,m),j=1,3),m=1,nres)
+ write (iout,'(8f10.5)') ((c(j,m+nres),j=1,3),m=nnt,nct)
write (iout,*) "Internal coordinates"
write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
#endif
do j=1,nss
+ if (dyn_ss) then
+ call xdrfint_(ixdrf, idssb(j), iret)
+ call xdrfint_(ixdrf, jdssb(j), iret)
+ idssb(j)=idssb(j)-nres
+ jdssb(j)=jdssb(j)-nres
+ else
call xdrfint_(ixdrf, ihpb(j), iret)
call xdrfint_(ixdrf, jhpb(j), iret)
+ endif
enddo
call xdrfint_(ixdrf, nprop, iret)
if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
#endif
do j=1,nss
+ if (dyn_ss) then
+ call xdrfint(ixdrf, idssb(j), iret)
+ call xdrfint(ixdrf, jdssb(j), iret)
+cc idssb(j)=idssb(j)-nres
+cc jdssb(j)=jdssb(j)-nres
+cc write(iout,*) idssb(j),jdssb(j)
+ else
call xdrfint(ixdrf, ihpb(j), iret)
call xdrfint(ixdrf, jhpb(j), iret)
+ endif
enddo
call xdrfint(ixdrf, nprop, iret)
c write (iout,*) "nprop",nprop
& ((csingle(l,k+nres),l=1,3),k=nnt,nct),
& nss,(ihpb(k),jhpb(k),k=1,nss),
& eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
+cc write(iout,*), 'NAWEJ',i,eini
if (indpdb.gt.0) then
do k=1,nres
do l=1,3
& " the value read in: ",energia(0),eini," point",
& iii+1,indstart(me1)+iii," T",
& 1.0d0/(1.987D-3*beta_h(ib,ipar))
+ call enerprint(energia(0),fT)
+ call pdbout(iii+1,beta_h(ib,ipar),
+ & eini,energia(0),0.0d0,rmsdev)
+ write (iout,*)
+
errmsg_count=errmsg_count+1
if (errmsg_count.gt.maxerrmsg_count)
& write (iout,*) "Too many warning messages"
do k=1,21
enetb(k,iii+1,iparm)=energia(k)
enddo
+c write(iout,*) "iCHUJ TU STRZELI",i,iii,entfac(i)
+c call enerprint(energia(0),fT)
#ifdef DEBUG
write (iout,'(2i5,f10.1,3e15.5)') i,iii,
& 1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
#else
do i=1,ntot(islice)
#endif
+cc if (dyn_ss) then
+cc read(ientout,rec=i,err=101)
+cc & ((csingle(l,k),l=1,3),k=1,nres),
+cc & ((csingle(l,k+nres),l=1,3),k=nnt,nct),
+cc & nss,(idssb(k),jdssb(k),k=1,nss),
+cc & eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
+cc idssb(k)=idssb(k)-nres
+cc jdssb(k)=jdssb(k)-nres
+cc else
read(ientout,rec=i,err=101)
& ((csingle(l,k),l=1,3),k=1,nres),
& ((csingle(l,k+nres),l=1,3),k=nnt,nct),
& nss,(ihpb(k),jhpb(k),k=1,nss),
& eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
-c write (iout,*) iR,ib,iparm,eini,efree
+cc endif
+cc write (iout,*) 'CC', iR,ib,iparm,eini,efree
do j=1,2*nres
do k=1,3
c(k,j)=csingle(k,j)
iscore=0
if (indpdb.gt.0) then
call conf_compar(i,.false.,.true.)
- endif
+ endif
+c if (dyn_ss) then
if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i)
& ((csingle(l,k),l=1,3),k=1,nres),
& ((csingle(l,k+nres),l=1,3),k=nnt,nct),
& nss,(ihpb(k),jhpb(k),k=1,nss),
c & potE(i,iparm),-entfac(i),rms_nat,iscore
& potE(i,nparmset),-entfac(i),rms_nat,iscore
-c write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
+c else
+ if (bxfile .or.cxfile .or. ensembles.gt.0) write
+ & (ientin,rec=i)
+ & ((csingle(l,k),l=1,3),k=1,nres),
+ & ((csingle(l,k+nres),l=1,3),k=nnt,nct),
+ & nss,(ihpb(k),jhpb(k),k=1,nss),
+c & potE(i,iparm),-entfac(i),rms_nat,iscore
+ & potE(i,nparmset),-entfac(i),rms_nat,iscore
+c endif
+ write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
#ifndef MPI
if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),
& -entfac(i),rms_nat,iscore)
c call flush(iout)
do i=indstart(j),indend(j)
iii = iii+1
+cc if (dyn_ss) then
+cc read(ientin,rec=iii,err=101)
+cc & ((csingle(l,k),l=1,3),k=1,nres),
+cc & ((csingle(l,k+nres),l=1,3),k=nnt,nct),
+cc & nss,(idssb(k),jdssb(k),k=1,nss),
+cc & eini,efree,rmsdev,iscor
+cc idssb(k)=idssb(k)-nres
+cc jdssb(k)=jdssb(k)-nres
+cc else
read(ientin,rec=iii,err=101)
& ((csingle(l,k),l=1,3),k=1,nres),
& ((csingle(l,k+nres),l=1,3),k=nnt,nct),
& nss,(ihpb(k),jhpb(k),k=1,nss),
& eini,efree,rmsdev,iscor
+cc endif
if (bxfile .or. ensembles.gt.0) then
- write (ientout,rec=i)
+cc if (dyn_ss) then
+cc write (ientout,rec=i)
+cc & ((csingle(l,k),l=1,3),k=1,nres),
+cc & ((csingle(l,k+nres),l=1,3),k=nnt,nct),
+cc & nss,(idssb(k)+nres,jdssb(k)+nres,k=1,nss),
+cc & eini,efree,rmsdev,iscor
+cc else
+ write (ientout,rec=i)
& ((csingle(l,k),l=1,3),k=1,nres),
& ((csingle(l,k+nres),l=1,3),k=nnt,nct),
& nss,(ihpb(k),jhpb(k),k=1,nss),
& eini,efree,rmsdev,iscor
+cc write(iout,*) "W poszukiwaniu zlotych galotow"
+cc write(iout,*) "efree=",efree,iii
+cc endif
endif
if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
#ifdef DEBUG
c write (iout,*) "itmp",itmp
c call flush(iout)
+c write (iout,*) "CNZ",eini,dyn_ss
#if (defined(AIX) && !defined(JUBL))
call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
c call flush(iout)
call xdrfint_(ixdrf, nss, iret)
do j=1,nss
+cc if (dyn_ss) then
+cc call xdrfint_(ixdrf, idssb(j)+nres, iret)
+cc call xdrfint_(ixdrf, jdssb(j)+nres, iret)
+cc else
call xdrfint_(ixdrf, ihpb(j), iret)
call xdrfint_(ixdrf, jhpb(j), iret)
+cc endif
enddo
call xdrffloat_(ixdrf,real(eini),iret)
call xdrffloat_(ixdrf,real(efree),iret)
call xdrfint(ixdrf, nss, iret)
do j=1,nss
+cc if (dyn_ss) then
+cc call xdrfint(ixdrf, idssb(j), iret)
+cc call xdrfint(ixdrf, jdssb(j), iret)
+cc idssb(j)=idssb(j)-nres
+cc jdssb(j)=jdssb(j)-nres
+cc else
call xdrfint(ixdrf, ihpb(j), iret)
call xdrfint(ixdrf, jhpb(j), iret)
+cc endif
enddo
call xdrffloat(ixdrf,real(eini),iret)
- call xdrffloat(ixdrf,real(efree),iret)
+ call xdrffloat(ixdrf,real(efree),iret)
call xdrffloat(ixdrf,real(rmsdev),iret)
call xdrfint(ixdrf,iscor,iret)
#endif
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
& +wvdwpp*evdw1
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
energia(19)=esccor
energia(20)=edihcnstr
energia(21)=evdw_t
+c if (dyn_ss) call dyn_set_nss
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
include 'COMMON.ENEPS'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SBRIDGE'
logical lprn
common /srutu/icall
integer icant
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
+C in case of diagnostics write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j)
+C /06/28/2013 Adasko: In case of dyn_ss - dynamic disulfide bond
+C formation no electrostatic interactions should be calculated. If it
+C would be allowed NaN would appear
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+C /06/28/2013 Adasko: dyn_ss_mask is logical statement wheather this Cys
+C residue can or cannot form disulfide bond. There is still bug allowing
+C Cys...Cys...Cys bond formation
+ call dyn_ssbond_ene(i,j,evdwij)
+C /06/28/2013 Adasko: dyn_ssbond_ene is dynamic SS bond foration energy
+C function in ssMD.F
+ evdw=evdw+evdwij
+c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+c & 'evdw',i,j,evdwij,' ss'
+ ELSE
ind=ind+1
itypj=itype(j)
dscj_inv=vbld_inv(j+nres)
c write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj,
c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
c & aux*e2/eps(itypi,itypj)
+c write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
C Calculate angular part of the gradient.
call sc_grad
endif
+ ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
+ if (.not.dyn_ss .and. i.le.nss) then
+C 15/02/13 CC dynamic SSbond - additional check
if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
+ endif
cd write (iout,*) "eij",eij
else if (ii.gt.nres .and. jj.gt.nres) then
c Restraints from contact prediction
deltat12=om2-om1+2.0d0
cosphi=om12-om1*om2
eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
- & +akct*deltad*deltat12
+ & +akct*deltad*deltat12+ebr
+c & +akct*deltad*deltat12
& +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
-c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
-c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
-c & " deltat12",deltat12," eij",eij
+ write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
+ & " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
+ & " deltat12",deltat12," eij",eij,"ebr",ebr
ed=2*akcm*deltad+akct*deltat12
pom1=akct*deltad
pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
include 'COMMON.FFIELD'
include 'COMMON.CONTROL'
double precision u(3),ud(3)
+ logical :: lprn=.false.
estr=0.0d0
do i=nnt+1,nct
diff = vbld(i)-vbldp0
nbi=nbondterm(iti)
if (nbi.eq.1) then
diff=vbld(i+nres)-vbldsc0(1,iti)
- write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+ if (lprn)
+ & write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
& AKSC(1,iti),AKSC(1,iti)*diff*diff
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
do j=1,3
usum=usum+uprod1
usumsqder=usumsqder+ud(j)*uprod2
enddo
- write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
+ if (lprn)
+ & write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
& AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
estr=estr+uprod/usum
do j=1,3
escloc = escloc + sumene
c write (2,*) "escloc",escloc
if (.not. calc_grad) goto 1
-#ifdef DEBUG
+
+#ifdef DEBUG2
C
C This section to check the numerical derivatives of the energy of ith side
C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
write (ipdb,30) ica(nct),ica(nct)+1
endif
do i=1,nss
+c if(dyn_ss) then
+c write (ipdb,30) ica(idssb(i))+1,ica(jdssb(i))+1
+c else
write (ipdb,30) ica(ihpb(i))+1,ica(jhpb(i))+1
+c endif
enddo
10 FORMAT ('ATOM',I7,' CA ',A3,I6,4X,3F8.3)
20 FORMAT ('ATOM',I7,' CB ',A3,I6,4X,3F8.3)
- double precision ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,dhpb,
- & dhpb1,forcon,weidis
- integer ns,nss,nfree,iss,ihpb,jhpb,nhpb,link_start,link_end,
- & ibecarb
- common /sbridge/ ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,ns,nss,
- & nfree,iss(maxss)
+ double precision ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss
+ integer ns,nss,nfree,iss
+ common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,
+ & ns,nss,nfree,iss(maxss)
+ double precision dhpb,dhpb1,forcon
+ integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb
common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
& ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb
+ double precision weidis
common /restraints/ weidis
+ integer link_start,link_end
common /links_split/ link_start,link_end
+ double precision Ht,dyn_ssbond_ij
+ logical dyn_ss,dyn_ss_mask
+ common /dyn_ssbond/ dyn_ssbond_ij(maxres,maxres),
+ & idssb(maxdim),jdssb(maxdim),
+ & Ht,dyn_ss,dyn_ss_mask(maxres)
--- /dev/null
+ integer king,idint,idreal,idchar,is_done
+ parameter (king=0,idint=1105,idreal=1729,idchar=1597,is_done=1)
+ integer me,cg_rank,fg_rank,fg_rank1,nodes,Nprocs,nfgtasks,kolor,
+ & koniec(0:maxprocs-1),WhatsUp,ifinish(maxprocs-1),CG_COMM,FG_COMM,
+ & FG_COMM1,CONT_FROM_COMM,CONT_TO_COMM,lentyp(0:maxprocs-1),
+ & kolor1,key1,nfgtasks1,MyRank,
+ & max_gs_size
+ logical yourjob, finished, cgdone
+ common/setup/me,MyRank,cg_rank,fg_rank,fg_rank1,nodes,Nprocs,
+ & nfgtasks,nfgtasks1,
+ & max_gs_size,kolor,koniec,WhatsUp,ifinish,CG_COMM,FG_COMM,
+ & FG_COMM1,CONT_FROM_COMM,CONT_TO_COMM,lentyp
+ integer MPI_UYZ,MPI_UYZGRAD,MPI_MU,MPI_MAT1,MPI_MAT2,
+ & MPI_THET,MPI_GAM,
+ & MPI_ROTAT1(0:1),MPI_ROTAT2(0:1),MPI_ROTAT_OLD(0:1),
+ & MPI_PRECOMP11(0:1),MPI_PRECOMP12(0:1),MPI_PRECOMP22(0:1),
+ & MPI_PRECOMP23(0:1)
+ common /types/ MPI_UYZ,MPI_UYZGRAD,MPI_MU,MPI_MAT1,MPI_MAT2,
+ & MPI_THET,MPI_GAM,
+ & MPI_ROTAT1,MPI_ROTAT2,MPI_ROTAT_OLD,MPI_PRECOMP11,MPI_PRECOMP12,
+ & MPI_PRECOMP22,MPI_PRECOMP23
ihpb(i)=0
jhpb(i)=0
enddo
+ do i=1,maxres
+ dyn_ss_mask(i)=.false.
+ enddo
C
C Initialize timing.
C
cd & (ihpb(i),jhpb(i),i=1,nss)
do i=nnt,nct-1
scheck=.false.
+ if (dyn_ss) go to 10
do ii=1,nss
if (ihpb(ii).eq.i+nres) then
scheck=.true.
endif
enddo ! i
#endif
- if (lprint) then
+ if (lprint) then
write (iout,'(a)') 'SC-p interaction array:'
do i=iatscp_s,iatscp_e
write (iout,'(i3,2(2x,2i3))')
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& beta_h(ib,iparm)*etot-entfac(i)
potE(i,iparm)=etot
#ifdef DEBUG
- write (iout,*) i,indstart(me)+i-1,ib,
+ write (iout,*) 'EEE',i,indstart(me)+i-1,ib,
& 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
& -entfac(i),Fdimless(i)
#endif
enddo
eini=fdimless(i)
call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
+cc if (temper.eq.300.0d0) then
+cc write(iout,*) 'Gral i=',iperm(i) ,eini,enepot(i),efree
+cc flush(iout)
+cc endif
enddo
+cc if (temper.eq.300.0d0) then
+cc write(iout,*) 'Gral i=',i ,eini,enepot(i),efree
+cc flush(iout)
+cc endif
#ifdef MPI
endif
#endif
enddo
endif
endif
+ if (ns.gt.0.and.dyn_ss) then
+C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
+C the bond
+ do i=nss+1,nhpb
+C /06/28/2013 Adasko: nss number of full SS bonds
+ ihpb(i-nss)=ihpb(i)
+ jhpb(i-nss)=jhpb(i)
+ forcon(i-nss)=forcon(i)
+ dhpb(i-nss)=dhpb(i)
+ enddo
+ nhpb=nhpb-nss
+ nss=0
+ call hpb_partition
+ do i=1,ns
+ dyn_ss_mask(iss(i))=.true.
+C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
+c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
+ enddo
+ endif
return
end
c------------------------------------------------------------------------------
key = wname(i)(:ilen(wname(i)))
call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
enddo
-
+ call reada(controlcard,"D0CM",d0cm,3.78d0)
+ call reada(controlcard,"AKCM",akcm,15.1d0)
+ call reada(controlcard,"AKTH",akth,11.0d0)
+ call reada(controlcard,"AKCT",akct,12.0d0)
+ call reada(controlcard,"V1SS",v1ss,-1.08d0)
+ call reada(controlcard,"V2SS",v2ss,7.61d0)
+ call reada(controlcard,"V3SS",v3ss,13.7d0)
+ call reada(controlcard,"EBR",ebr,-5.50D0)
+c dyn_ss=(index(controlcard,'DYN_SS').gt.0)
write (iout,*) "iparm",iparm," myparm",myparm
-c If reading not own parameters, skip assignment
+c do i=1,maxres
+c dyn_ss_mask(i)=.false.
+c enddo
+ do i=1,maxres-1
+ do j=i+1,maxres
+ dyn_ssbond_ij(i,j)=1.0d300
+ enddo
+ enddo
+ call reada(controlcard,"HT",Ht,0.0D0)
+c if(me.eq.king.or..not.out1file) then
+c print *,'indpdb=',indpdb,' pdbref=',pdbref
+c endif
+c If reading not own parameters, skip assignment
+cc write(iout,*) "KURWA", ww(15)
if (iparm.eq.myparm .or. .not.separate_parset) then
c
wsccor=ww(19)
endif
+cc write(iout,*) "KURWA", wstrain,akcm,akth,wsc,dyn_ss
call card_concat(controlcard,.false.)
enddo
enddo
C
-C Define the SC-p interaction constants
-C
+C Define the SC-p interaction constants and SS bond potentials
+C
+ if (dyn_ss) then
+ ss_depth=ebr/wsc-0.25*eps(1,1)
+ Ht=Ht/wsc-0.25*eps(1,1)
+ akcm=akcm*wstrain/wsc
+ akth=akth*wstrain/wsc
+ akct=akct*wstrain/wsc
+ v1ss=v1ss*wstrain/wsc
+ v2ss=v2ss*wstrain/wsc
+ v3ss=v3ss*wstrain/wsc
+ else
+ ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+ endif
+ write (iout,*) "Parameters of the SS-bond potential:"
+ write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
+ & " AKCT",akct
+ write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
+ write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
+ write (iout,*)" HT",Ht
+
#ifdef OLDSCP
do i=1,20
C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
C
C Define the constants of the disulfide bridge
C
- ebr=-5.50D0
+c ebr=-5.50D0
c
c Old arbitrary potential - commented out.
c
c energy surface of diethyl disulfide.
c A. Liwo and U. Kozlowska, 11/24/03
c
- D0CM = 3.78d0
- AKCM = 15.1d0
- AKTH = 11.0d0
- AKCT = 12.0d0
- V1SS =-1.08d0
- V2SS = 7.61d0
- V3SS = 13.7d0
+c D0CM = 3.78d0
+c AKCM = 15.1d0
+c AKTH = 11.0d0
+c AKCT = 12.0d0
+c V1SS =-1.08d0
+c V2SS = 7.61d0
+c V3SS = 13.7d0
- if (lprint) then
- write (iout,'(/a)') "Disulfide bridge parameters:"
- write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
- write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
- write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
- write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
- & ' v3ss:',v3ss
- endif
+c if (lprint) then
+c write (iout,'(/a)') "Disulfide bridge parameters:"
+c write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
+c write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
+c write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
+c write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
+c & ' v3ss:',v3ss
+c endif
return
end
include "COMMON.FREE"
include "COMMON.CONTROL"
include "COMMON.ENERGIES"
+ include "COMMON.SBRIDGE"
character*800 controlcard
integer i,j,k,ii,n_ene_found
integer ind,itype1,itype2,itypf,itypsc,itypp
& " CONSTR_DIST",constr_dist
refstr = index(controlcard,'REFSTR').gt.0
pdbref = index(controlcard,'PDBREF').gt.0
+ dyn_ss=(index(controlcard,'DYN_SS').gt.0)
+C /06/28/2013 Adasko: dyn_ss is keyword allowing to break and create bond
+C disulfide bond. Note that in conterary to dynamics this in
+C CONTROLCARD. The bond are read in molread_zs.F
call flush(iout)
return
end
--- /dev/null
+c----------------------------------------------------------------------------
+ subroutine check_energies
+c implicit none
+
+c Includes
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.LOCAL'
+ include 'COMMON.GEO'
+
+c External functions
+ double precision ran_number
+ external ran_number
+
+c Local variables
+ integer i,j,k,l,lmax,p,pmax
+ double precision rmin,rmax
+ double precision eij
+
+ double precision d
+ double precision wi,rij,tj,pj
+
+
+c return
+
+ i=5
+ j=14
+
+ d=dsc(1)
+ rmin=2.0D0
+ rmax=12.0D0
+
+ lmax=10000
+ pmax=1
+
+ do k=1,3
+ c(k,i)=0.0D0
+ c(k,j)=0.0D0
+ c(k,nres+i)=0.0D0
+ c(k,nres+j)=0.0D0
+ enddo
+
+ do l=1,lmax
+
+ct wi=ran_number(0.0D0,pi)
+c wi=ran_number(0.0D0,pi/6.0D0)
+c wi=0.0D0
+ct tj=ran_number(0.0D0,pi)
+ct pj=ran_number(0.0D0,pi)
+c pj=ran_number(0.0D0,pi/6.0D0)
+c pj=0.0D0
+
+ do p=1,pmax
+ct rij=ran_number(rmin,rmax)
+
+ c(1,j)=d*sin(pj)*cos(tj)
+ c(2,j)=d*sin(pj)*sin(tj)
+ c(3,j)=d*cos(pj)
+
+ c(3,nres+i)=-rij
+
+ c(1,i)=d*sin(wi)
+ c(3,i)=-rij-d*cos(wi)
+
+ do k=1,3
+ dc(k,nres+i)=c(k,nres+i)-c(k,i)
+ dc_norm(k,nres+i)=dc(k,nres+i)/d
+ dc(k,nres+j)=c(k,nres+j)-c(k,j)
+ dc_norm(k,nres+j)=dc(k,nres+j)/d
+ enddo
+
+ call dyn_ssbond_ene(i,j,eij)
+ enddo
+ enddo
+
+ call exit(1)
+
+ return
+ end
+
+C-----------------------------------------------------------------------------
+
+ subroutine dyn_ssbond_ene(resi,resj,eij)
+c implicit none
+
+c Includes
+ include 'DIMENSIONS'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+ include 'COMMON.MD'
+#endif
+#endif
+
+c External functions
+ double precision h_base
+ external h_base
+
+c Input arguments
+ integer resi,resj
+
+c Output arguments
+ double precision eij
+
+c Local variables
+ logical havebond
+c integer itypi,itypj,k,l
+ double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
+ double precision sig0ij,ljd,sig,fac,e1,e2
+ double precision dcosom1(3),dcosom2(3),ed
+ double precision pom1,pom2
+ double precision ljA,ljB,ljXs
+ double precision d_ljB(1:3)
+ double precision ssA,ssB,ssC,ssXs
+ double precision ssxm,ljxm,ssm,ljm
+ double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+ double precision f1,f2,h1,h2,hd1,hd2
+ double precision omega,delta_inv,deltasq_inv,fac1,fac2
+c-------FIRST METHOD
+ double precision xm,d_xm(1:3)
+c-------END FIRST METHOD
+c-------SECOND METHOD
+c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
+c-------END SECOND METHOD
+
+c-------TESTING CODE
+ logical checkstop,transgrad
+ common /sschecks/ checkstop,transgrad
+
+ integer icheck,nicheck,jcheck,njcheck
+ double precision echeck(-1:1),deps,ssx0,ljx0
+c-------END TESTING CODE
+
+
+ i=resi
+ j=resj
+
+ itypi=itype(i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+
+ itypj=itype(j)
+ xj=c(1,nres+j)-c(1,nres+i)
+ yj=c(2,nres+j)-c(2,nres+i)
+ zj=c(3,nres+j)-c(3,nres+i)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ dscj_inv=vbld_inv(j+nres)
+
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
+c The following are set in sc_angular
+c erij(1)=xj*rij
+c erij(2)=yj*rij
+c erij(3)=zj*rij
+c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+c om12=dxi*dxj+dyi*dyj+dzi*dzj
+ call sc_angular
+ rij=1.0D0/rij ! Reset this so it makes sense
+
+ sig0ij=sigma(itypi,itypj)
+ sig=sig0ij*dsqrt(1.0D0/sigsq)
+
+ ljXs=sig-sig0ij
+ ljA=eps1*eps2rt**2*eps3rt**2
+ ljB=ljA*bb(itypi,itypj)
+ ljA=ljA*aa(itypi,itypj)
+ ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+
+ ssXs=d0cm
+ deltat1=1.0d0-om1
+ deltat2=1.0d0+om2
+ deltat12=om2-om1+2.0d0
+ cosphi=om12-om1*om2
+ ssA=akcm
+ ssB=akct*deltat12
+ ssC=ss_depth
+ & +akth*(deltat1*deltat1+deltat2*deltat2)
+ & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
+ ssxm=ssXs-0.5D0*ssB/ssA
+
+c-------TESTING CODE
+c$$$c Some extra output
+c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
+c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
+c$$$ if (ssx0.gt.0.0d0) then
+c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
+c$$$ else
+c$$$ ssx0=ssxm
+c$$$ endif
+c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
+c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
+c$$$ return
+c-------END TESTING CODE
+
+c-------TESTING CODE
+c Stop and plot energy and derivative as a function of distance
+ if (checkstop) then
+ ssm=ssC-0.25D0*ssB*ssB/ssA
+ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+ if (ssm.lt.ljm .and.
+ & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
+ nicheck=1000
+ njcheck=1
+ deps=0.5d-7
+ else
+ checkstop=.false.
+ endif
+ endif
+ if (.not.checkstop) then
+ nicheck=0
+ njcheck=-1
+ endif
+
+ do icheck=0,nicheck
+ do jcheck=-1,njcheck
+ if (checkstop) rij=(ssxm-1.0d0)+
+ & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
+c-------END TESTING CODE
+
+ if (rij.gt.ljxm) then
+ havebond=.false.
+ ljd=rij-ljXs
+ fac=(1.0D0/ljd)**expon
+ e1=fac*fac*aa(itypi,itypj)
+ e2=fac*bb(itypi,itypj)
+ eij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=eij*eps3rt
+ eps3der=eij*eps2rt
+ eij=eij*eps2rt*eps3rt
+
+ sigder=-sig/sigsq
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ ed=-expon*(e1+eij)/ljd
+ sigder=ed*sigder
+ eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
+ eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
+ eom12=eij*eps1_om12+eps2der*eps2rt_om12
+ & -2.0D0*alf12*eps3der+sigder*sigsq_om12
+ else if (rij.lt.ssxm) then
+ havebond=.true.
+ ssd=rij-ssXs
+ eij=ssA*ssd*ssd+ssB*ssd+ssC
+
+ ed=2*akcm*ssd+akct*deltat12
+ pom1=akct*ssd
+ pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
+ eom1=-2*akth*deltat1-pom1-om2*pom2
+ eom2= 2*akth*deltat2+pom1-om1*pom2
+ eom12=pom2
+ else
+ omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
+
+ d_ssxm(1)=0.5D0*akct/ssA
+ d_ssxm(2)=-d_ssxm(1)
+ d_ssxm(3)=0.0D0
+
+ d_ljxm(1)=sig0ij/sqrt(sigsq**3)
+ d_ljxm(2)=d_ljxm(1)*sigsq_om2
+ d_ljxm(3)=d_ljxm(1)*sigsq_om12
+ d_ljxm(1)=d_ljxm(1)*sigsq_om1
+
+c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
+ xm=0.5d0*(ssxm+ljxm)
+ do k=1,3
+ d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
+ enddo
+ if (rij.lt.xm) then
+ havebond=.true.
+ ssm=ssC-0.25D0*ssB*ssB/ssA
+ d_ssm(1)=0.5D0*akct*ssB/ssA
+ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
+ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
+ d_ssm(3)=omega
+ f1=(rij-xm)/(ssxm-xm)
+ f2=(rij-ssxm)/(xm-ssxm)
+ h1=h_base(f1,hd1)
+ h2=h_base(f2,hd2)
+ eij=ssm*h1+Ht*h2
+ delta_inv=1.0d0/(xm-ssxm)
+ deltasq_inv=delta_inv*delta_inv
+ fac=ssm*hd1-Ht*hd2
+ fac1=deltasq_inv*fac*(xm-rij)
+ fac2=deltasq_inv*fac*(rij-ssxm)
+ ed=delta_inv*(Ht*hd2-ssm*hd1)
+ eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
+ eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
+ eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
+ else
+ havebond=.false.
+ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+ d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
+ d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
+ d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
+ + alf12/eps3rt)
+ d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
+ f1=(rij-ljxm)/(xm-ljxm)
+ f2=(rij-xm)/(ljxm-xm)
+ h1=h_base(f1,hd1)
+ h2=h_base(f2,hd2)
+ eij=Ht*h1+ljm*h2
+ delta_inv=1.0d0/(ljxm-xm)
+ deltasq_inv=delta_inv*delta_inv
+ fac=Ht*hd1-ljm*hd2
+ fac1=deltasq_inv*fac*(ljxm-rij)
+ fac2=deltasq_inv*fac*(rij-xm)
+ ed=delta_inv*(ljm*hd2-Ht*hd1)
+ eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
+ eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
+ eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
+ endif
+c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
+
+c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
+c$$$ ssd=rij-ssXs
+c$$$ ljd=rij-ljXs
+c$$$ fac1=rij-ljxm
+c$$$ fac2=rij-ssxm
+c$$$
+c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
+c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
+c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
+c$$$
+c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
+c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
+c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
+c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
+c$$$ d_ssm(3)=omega
+c$$$
+c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
+c$$$ do k=1,3
+c$$$ d_ljm(k)=ljm*d_ljB(k)
+c$$$ enddo
+c$$$ ljm=ljm*ljB
+c$$$
+c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
+c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
+c$$$ d_ss(2)=akct*ssd
+c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
+c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
+c$$$ d_ss(3)=omega
+c$$$
+c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
+c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
+c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
+c$$$ do k=1,3
+c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
+c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
+c$$$ enddo
+c$$$ ljf=ljm+ljf*ljB*fac1*fac1
+c$$$
+c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
+c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
+c$$$ h1=h_base(f1,hd1)
+c$$$ h2=h_base(f2,hd2)
+c$$$ eij=ss*h1+ljf*h2
+c$$$ delta_inv=1.0d0/(ljxm-ssxm)
+c$$$ deltasq_inv=delta_inv*delta_inv
+c$$$ fac=ljf*hd2-ss*hd1
+c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
+c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
+c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
+c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
+c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
+c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
+c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
+c$$$
+c$$$ havebond=.false.
+c$$$ if (ed.gt.0.0d0) havebond=.true.
+c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
+
+ endif
+
+ if (havebond) then
+#ifndef CLUST
+#ifndef WHAM
+c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
+c write(iout,'(a15,f12.2,f8.1,2i5)')
+c & "SSBOND_E_FORM",totT,t_bath,i,j
+c endif
+#endif
+#endif
+ write(iout,*), 'DYN_SS_BOND',i,j,eij
+ dyn_ssbond_ij(i,j)=eij
+ else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
+ dyn_ssbond_ij(i,j)=1.0d300
+#ifndef CLUST
+#ifndef WHAM
+c write(iout,'(a15,f12.2,f8.1,2i5)')
+c & "SSBOND_E_BREAK",totT,t_bath,i,j
+#endif
+#endif
+ endif
+
+c-------TESTING CODE
+ if (checkstop) then
+ if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
+ & "CHECKSTOP",rij,eij,ed
+ echeck(jcheck)=eij
+ endif
+ enddo
+ if (checkstop) then
+ write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
+ endif
+ enddo
+ if (checkstop) then
+ transgrad=.true.
+ checkstop=.false.
+ endif
+c-------END TESTING CODE
+
+ do k=1,3
+ dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
+ dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
+ enddo
+ do k=1,3
+ gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ enddo
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
+ & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
+ & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+cgrad do k=i,j-1
+cgrad do l=1,3
+cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
+cgrad enddo
+cgrad enddo
+
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+
+ return
+ end
+
+C-----------------------------------------------------------------------------
+
+ double precision function h_base(x,deriv)
+c A smooth function going 0->1 in range [0,1]
+c It should NOT be called outside range [0,1], it will not work there.
+ implicit none
+
+c Input arguments
+ double precision x
+
+c Output arguments
+ double precision deriv
+
+c Local variables
+ double precision xsq
+
+
+c Two parabolas put together. First derivative zero at extrema
+c$$$ if (x.lt.0.5D0) then
+c$$$ h_base=2.0D0*x*x
+c$$$ deriv=4.0D0*x
+c$$$ else
+c$$$ deriv=1.0D0-x
+c$$$ h_base=1.0D0-2.0D0*deriv*deriv
+c$$$ deriv=4.0D0*deriv
+c$$$ endif
+
+c Third degree polynomial. First derivative zero at extrema
+ h_base=x*x*(3.0d0-2.0d0*x)
+ deriv=6.0d0*x*(1.0d0-x)
+
+c Fifth degree polynomial. First and second derivatives zero at extrema
+c$$$ xsq=x*x
+c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
+c$$$ deriv=x-1.0d0
+c$$$ deriv=deriv*deriv
+c$$$ deriv=30.0d0*xsq*deriv
+
+ return
+ end
+
+c----------------------------------------------------------------------------
+
+ subroutine dyn_set_nss
+c Adjust nss and other relevant variables based on dyn_ssbond_ij
+c implicit none
+
+c Includes
+ include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+#endif
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+c#ifdef CLUST
+c#ifdef WHAM
+ integer max_fg_procs,maxprocs
+ parameter (max_fg_procs=128)
+ parameter (maxprocs=128)
+c#endif
+c#endif
+cc include 'COMMON.SETUP'
+#ifndef CLUST
+#ifndef WHAM
+ include 'COMMON.MD'
+ include 'COMMON.SETUP'
+#endif
+#endif
+
+c Local variables
+ double precision emin
+ integer i,j,imin
+ integer diff,allflag(maxdim),allnss,
+ & allihpb(maxdim),alljhpb(maxdim),
+ & newnss,newihpb(maxdim),newjhpb(maxdim)
+ logical found
+ integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
+ integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
+
+ allnss=0
+ do i=1,nres-1
+ do j=i+1,nres
+ if (dyn_ssbond_ij(i,j).lt.1.0d300) then
+ allnss=allnss+1
+ allflag(allnss)=0
+ allihpb(allnss)=i
+ alljhpb(allnss)=j
+ endif
+ enddo
+ enddo
+
+cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
+
+ 1 emin=1.0d300
+ do i=1,allnss
+ if (allflag(i).eq.0 .and.
+ & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
+ emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
+ imin=i
+ endif
+ enddo
+ if (emin.lt.1.0d300) then
+ allflag(imin)=1
+ do i=1,allnss
+ if (allflag(i).eq.0 .and.
+ & (allihpb(i).eq.allihpb(imin) .or.
+ & alljhpb(i).eq.allihpb(imin) .or.
+ & allihpb(i).eq.alljhpb(imin) .or.
+ & alljhpb(i).eq.alljhpb(imin))) then
+ allflag(i)=-1
+ endif
+ enddo
+ goto 1
+ endif
+
+cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
+
+ newnss=0
+ do i=1,allnss
+ if (allflag(i).eq.1) then
+ newnss=newnss+1
+ newihpb(newnss)=allihpb(i)
+ newjhpb(newnss)=alljhpb(i)
+ endif
+ enddo
+
+#ifdef MPI
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(newnss,g_newnss,1,
+ & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+ call MPI_Gather(newnss,1,MPI_INTEGER,
+ & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_newnss(i-1)+displ(i-1)
+ enddo
+ call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
+ & g_newihpb,i_newnss,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
+ & g_newjhpb,i_newnss,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ if(fg_rank.eq.0) then
+c print *,'g_newnss',g_newnss
+c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
+c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
+ newnss=g_newnss
+ do i=1,newnss
+ newihpb(i)=g_newihpb(i)
+ newjhpb(i)=g_newjhpb(i)
+ enddo
+ endif
+ endif
+#endif
+
+ diff=newnss-nss
+
+cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
+
+ do i=1,nss
+ found=.false.
+ do j=1,newnss
+ if (idssb(i).eq.newihpb(j) .and.
+ & jdssb(i).eq.newjhpb(j)) found=.true.
+ enddo
+#ifndef CLUST
+#ifndef WHAM
+ if (.not.found.and.fg_rank.eq.0)
+ & write(iout,'(a15,f12.2,f8.1,2i5)')
+ & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
+#endif
+#endif
+ enddo
+
+ do i=1,newnss
+ found=.false.
+ do j=1,nss
+ if (newihpb(i).eq.idssb(j) .and.
+ & newjhpb(i).eq.jdssb(j)) found=.true.
+ enddo
+#ifndef CLUST
+#ifndef WHAM
+ if (.not.found.and.fg_rank.eq.0)
+ & write(iout,'(a15,f12.2,f8.1,2i5)')
+ & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
+#endif
+#endif
+ enddo
+
+ nss=newnss
+ do i=1,nss
+ idssb(i)=newihpb(i)
+ jdssb(i)=newjhpb(i)
+ enddo
+
+ return
+ end
+
+c----------------------------------------------------------------------------
+#ifdef NIEWIADOM
+c#ifdef WHAM
+ subroutine read_ssHist
+ implicit none
+
+c Includes
+ include 'DIMENSIONS'
+ include "DIMENSIONS.FREE"
+ include 'COMMON.FREE'
+
+c Local variables
+ integer i,j
+ character*80 controlcard
+
+ do i=1,dyn_nssHist
+ call card_concat(controlcard,.true.)
+ read(controlcard,*)
+ & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
+ enddo
+
+ return
+ end
+#endif
+
+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-----------------------------------------------------------------------------
enddo
enddo
c Store disulfide-bond parameters
+ ht_all(iparm)=ht
+ ss_depth_all(iparm)=ss_depth
ebr_all(iparm)=ebr
d0cm_all(iparm)=d0cm
akcm_all(iparm)=akcm
enddo
enddo
c Restore disulfide-bond parameters
+ ht=ht_all(iparm)
+ ss_depth=ss_depth_all(iparm)
ebr=ebr_all(iparm)
d0cm=d0cm_all(iparm)
akcm=akcm_all(iparm)
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr