1 c----------------------------------------------------------------------------
2 subroutine check_energies
9 include 'COMMON.IOUNITS'
10 include 'COMMON.SBRIDGE'
11 include 'COMMON.LOCAL'
15 double precision ran_number
19 integer i,j,k,l,lmax,p,pmax
20 double precision rmin,rmax
24 double precision wi,rij,tj,pj
48 ct wi=ran_number(0.0D0,pi)
49 c wi=ran_number(0.0D0,pi/6.0D0)
51 ct tj=ran_number(0.0D0,pi)
52 ct pj=ran_number(0.0D0,pi)
53 c pj=ran_number(0.0D0,pi/6.0D0)
57 ct rij=ran_number(rmin,rmax)
59 c(1,j)=d*sin(pj)*cos(tj)
60 c(2,j)=d*sin(pj)*sin(tj)
69 dc(k,nres+i)=c(k,nres+i)-c(k,i)
70 dc_norm(k,nres+i)=dc(k,nres+i)/d
71 dc(k,nres+j)=c(k,nres+j)-c(k,j)
72 dc_norm(k,nres+j)=dc(k,nres+j)/d
75 call dyn_ssbond_ene(i,j,eij)
84 C-----------------------------------------------------------------------------
86 subroutine dyn_ssbond_ene(resi,resj,eij)
91 include 'COMMON.SBRIDGE'
92 include 'COMMON.CHAIN'
93 include 'COMMON.DERIV'
94 include 'COMMON.LOCAL'
95 include 'COMMON.INTERACT'
97 include 'COMMON.IOUNITS'
106 double precision h_base
117 c integer itypi,itypj,k,l
118 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
119 double precision sig0ij,ljd,sig,fac,e1,e2
120 double precision dcosom1(3),dcosom2(3),ed
121 double precision pom1,pom2
122 double precision ljA,ljB,ljXs
123 double precision d_ljB(1:3)
124 double precision ssA,ssB,ssC,ssXs
125 double precision ssxm,ljxm,ssm,ljm
126 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
127 double precision f1,f2,h1,h2,hd1,hd2
128 double precision omega,delta_inv,deltasq_inv,fac1,fac2
130 double precision xm,d_xm(1:3)
131 c-------END FIRST METHOD
132 c-------SECOND METHOD
133 c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
134 c-------END SECOND METHOD
137 logical checkstop,transgrad
138 common /sschecks/ checkstop,transgrad
140 integer icheck,nicheck,jcheck,njcheck
141 double precision echeck(-1:1),deps,ssx0,ljx0
142 c-------END TESTING CODE
149 dxi=dc_norm(1,nres+i)
150 dyi=dc_norm(2,nres+i)
151 dzi=dc_norm(3,nres+i)
152 dsci_inv=vbld_inv(i+nres)
155 xj=c(1,nres+j)-c(1,nres+i)
156 yj=c(2,nres+j)-c(2,nres+i)
157 zj=c(3,nres+j)-c(3,nres+i)
158 dxj=dc_norm(1,nres+j)
159 dyj=dc_norm(2,nres+j)
160 dzj=dc_norm(3,nres+j)
161 dscj_inv=vbld_inv(j+nres)
163 chi1=chi(itypi,itypj)
164 chi2=chi(itypj,itypi)
171 alf12=0.5D0*(alf1+alf2)
173 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
174 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
175 c The following are set in sc_angular
179 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
180 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
181 c om12=dxi*dxj+dyi*dyj+dzi*dzj
183 rij=1.0D0/rij ! Reset this so it makes sense
185 sig0ij=sigma(itypi,itypj)
186 sig=sig0ij*dsqrt(1.0D0/sigsq)
189 ljA=eps1*eps2rt**2*eps3rt**2
190 ljB=ljA*bb(itypi,itypj)
191 ljA=ljA*aa(itypi,itypj)
192 ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
197 deltat12=om2-om1+2.0d0
202 & +akth*(deltat1*deltat1+deltat2*deltat2)
203 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
204 ssxm=ssXs-0.5D0*ssB/ssA
207 c$$$c Some extra output
208 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
209 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
210 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
211 c$$$ if (ssx0.gt.0.0d0) then
212 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
216 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
217 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
218 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
220 c-------END TESTING CODE
223 c Stop and plot energy and derivative as a function of distance
225 ssm=ssC-0.25D0*ssB*ssB/ssA
226 ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
228 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
236 if (.not.checkstop) then
243 if (checkstop) rij=(ssxm-1.0d0)+
244 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
245 c-------END TESTING CODE
247 if (rij.gt.ljxm) then
250 fac=(1.0D0/ljd)**expon
251 e1=fac*fac*aa(itypi,itypj)
252 e2=fac*bb(itypi,itypj)
253 eij=eps1*eps2rt*eps3rt*(e1+e2)
256 eij=eij*eps2rt*eps3rt
259 e1=e1*eps1*eps2rt**2*eps3rt**2
260 ed=-expon*(e1+eij)/ljd
262 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
263 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
264 eom12=eij*eps1_om12+eps2der*eps2rt_om12
265 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
266 else if (rij.lt.ssxm) then
269 eij=ssA*ssd*ssd+ssB*ssd+ssC
271 ed=2*akcm*ssd+akct*deltat12
273 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
274 eom1=-2*akth*deltat1-pom1-om2*pom2
275 eom2= 2*akth*deltat2+pom1-om1*pom2
278 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
280 d_ssxm(1)=0.5D0*akct/ssA
284 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
285 d_ljxm(2)=d_ljxm(1)*sigsq_om2
286 d_ljxm(3)=d_ljxm(1)*sigsq_om12
287 d_ljxm(1)=d_ljxm(1)*sigsq_om1
289 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
292 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
296 ssm=ssC-0.25D0*ssB*ssB/ssA
297 d_ssm(1)=0.5D0*akct*ssB/ssA
298 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
299 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
301 f1=(rij-xm)/(ssxm-xm)
302 f2=(rij-ssxm)/(xm-ssxm)
306 delta_inv=1.0d0/(xm-ssxm)
307 deltasq_inv=delta_inv*delta_inv
309 fac1=deltasq_inv*fac*(xm-rij)
310 fac2=deltasq_inv*fac*(rij-ssxm)
311 ed=delta_inv*(Ht*hd2-ssm*hd1)
312 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
313 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
314 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
317 ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
318 d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
319 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
320 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
322 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
323 f1=(rij-ljxm)/(xm-ljxm)
324 f2=(rij-xm)/(ljxm-xm)
328 delta_inv=1.0d0/(ljxm-xm)
329 deltasq_inv=delta_inv*delta_inv
331 fac1=deltasq_inv*fac*(ljxm-rij)
332 fac2=deltasq_inv*fac*(rij-xm)
333 ed=delta_inv*(ljm*hd2-Ht*hd1)
334 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
335 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
336 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
338 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
340 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
346 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
347 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
348 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
350 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
351 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
352 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
353 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
356 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
358 c$$$ d_ljm(k)=ljm*d_ljB(k)
362 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
363 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
364 c$$$ d_ss(2)=akct*ssd
365 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
366 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
369 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
370 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
371 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
373 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
374 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
376 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
378 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
379 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
380 c$$$ h1=h_base(f1,hd1)
381 c$$$ h2=h_base(f2,hd2)
382 c$$$ eij=ss*h1+ljf*h2
383 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
384 c$$$ deltasq_inv=delta_inv*delta_inv
385 c$$$ fac=ljf*hd2-ss*hd1
386 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
387 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
388 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
389 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
390 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
391 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
392 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
394 c$$$ havebond=.false.
395 c$$$ if (ed.gt.0.0d0) havebond=.true.
396 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
403 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
404 c write(iout,'(a15,f12.2,f8.1,2i5)')
405 c & "SSBOND_E_FORM",totT,t_bath,i,j
409 dyn_ssbond_ij(i,j)=eij
410 else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
411 dyn_ssbond_ij(i,j)=1.0d300
414 c write(iout,'(a15,f12.2,f8.1,2i5)')
415 c & "SSBOND_E_BREAK",totT,t_bath,i,j
422 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
423 & "CHECKSTOP",rij,eij,ed
428 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
435 c-------END TESTING CODE
438 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
439 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
442 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
445 gvdwx(k,i)=gvdwx(k,i)-gg(k)
446 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
447 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
448 gvdwx(k,j)=gvdwx(k,j)+gg(k)
449 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
450 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
454 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
459 gvdwc(l,i)=gvdwc(l,i)-gg(l)
460 gvdwc(l,j)=gvdwc(l,j)+gg(l)
466 C-----------------------------------------------------------------------------
468 double precision function h_base(x,deriv)
469 c A smooth function going 0->1 in range [0,1]
470 c It should NOT be called outside range [0,1], it will not work there.
477 double precision deriv
483 c Two parabolas put together. First derivative zero at extrema
484 c$$$ if (x.lt.0.5D0) then
485 c$$$ h_base=2.0D0*x*x
489 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
490 c$$$ deriv=4.0D0*deriv
493 c Third degree polynomial. First derivative zero at extrema
494 h_base=x*x*(3.0d0-2.0d0*x)
495 deriv=6.0d0*x*(1.0d0-x)
497 c Fifth degree polynomial. First and second derivatives zero at extrema
499 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
501 c$$$ deriv=deriv*deriv
502 c$$$ deriv=30.0d0*xsq*deriv
507 c----------------------------------------------------------------------------
509 subroutine dyn_set_nss
510 c Adjust nss and other relevant variables based on dyn_ssbond_ij
515 include 'COMMON.SBRIDGE'
516 include 'COMMON.CHAIN'
517 include 'COMMON.IOUNITS'
525 double precision emin
527 integer diff,allflag(maxdim),allnss,
528 & allihpb(maxdim),alljhpb(maxdim),
529 & newnss,newihpb(maxdim),newjhpb(maxdim)
536 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
545 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
549 if (allflag(i).eq.0 .and.
550 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
551 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
555 if (emin.lt.1.0d300) then
558 if (allflag(i).eq.0 .and.
559 & (allihpb(i).eq.allihpb(imin) .or.
560 & alljhpb(i).eq.allihpb(imin) .or.
561 & allihpb(i).eq.alljhpb(imin) .or.
562 & alljhpb(i).eq.alljhpb(imin))) then
569 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
573 if (allflag(i).eq.1) then
575 newihpb(newnss)=allihpb(i)
576 newjhpb(newnss)=alljhpb(i)
581 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
586 if (ihpb(i).eq.newihpb(j) .and.
587 & jhpb(i).eq.newjhpb(j)) found=.true.
591 if (.not.found) write(iout,'(a15,f12.2,f8.1,2i5)')
592 & "SSBOND_BREAK",totT,t_bath,ihpb(i),jhpb(i)
600 if (newihpb(i).eq.ihpb(j) .and.
601 & newjhpb(i).eq.jhpb(j)) found=.true.
605 if (.not.found) write(iout,'(a15,f12.2,f8.1,2i5)')
606 & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
613 ihpb(nhpb+i)=ihpb(nss+i)
614 jhpb(nhpb+i)=jhpb(nss+i)
615 forcon(nhpb+i)=forcon(nss+i)
616 dhpb(nhpb+i)=dhpb(nss+i)
618 else if (diff.lt.0) then
620 ihpb(nss+i)=ihpb(nhpb+i)
621 jhpb(nss+i)=jhpb(nhpb+i)
622 forcon(nss+i)=forcon(nhpb+i)
623 dhpb(nss+i)=dhpb(nhpb+i)
637 c----------------------------------------------------------------------------
640 subroutine read_ssHist
645 include "DIMENSIONS.FREE"
646 include 'COMMON.FREE'
650 character*80 controlcard
653 call card_concat(controlcard,.true.)
655 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
662 c----------------------------------------------------------------------------
665 C-----------------------------------------------------------------------------
666 C-----------------------------------------------------------------------------
667 C-----------------------------------------------------------------------------
668 C-----------------------------------------------------------------------------
669 C-----------------------------------------------------------------------------
670 C-----------------------------------------------------------------------------
671 C-----------------------------------------------------------------------------
673 c$$$c-----------------------------------------------------------------------------
675 c$$$ subroutine ss_relax(i_in,j_in)
679 c$$$ include 'DIMENSIONS'
680 c$$$ include 'COMMON.VAR'
681 c$$$ include 'COMMON.CHAIN'
682 c$$$ include 'COMMON.IOUNITS'
683 c$$$ include 'COMMON.INTERACT'
685 c$$$c Input arguments
686 c$$$ integer i_in,j_in
688 c$$$c Local variables
689 c$$$ integer i,iretcode,nfun_sc
691 c$$$ double precision var(maxvar),e_sc,etot
698 c$$$ mask_side(i_in)=1
699 c$$$ mask_side(j_in)=1
701 c$$$c Minimize the two selected side-chains
702 c$$$ call overlap_sc(scfail) ! Better not fail!
703 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
710 c$$$c-------------------------------------------------------------
712 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
713 c$$$c Minimize side-chains only, starting from geom but without modifying
715 c$$$c If mask_r is already set, only the selected side-chains are minimized,
716 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
720 c$$$ include 'DIMENSIONS'
721 c$$$ include 'COMMON.IOUNITS'
722 c$$$ include 'COMMON.VAR'
723 c$$$ include 'COMMON.CHAIN'
724 c$$$ include 'COMMON.GEO'
725 c$$$ include 'COMMON.MINIM'
727 c$$$ common /srutu/ icall
729 c$$$c Output arguments
730 c$$$ double precision etot_sc
731 c$$$ integer iretcode,nfun
733 c$$$c External functions/subroutines
734 c$$$ external func_sc,grad_sc,fdum
736 c$$$c Local variables
738 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
740 c$$$ double precision rdum(1)
741 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
743 c$$$ integer i,nvar_restr
746 c$$$cmc start_minim=.true.
747 c$$$ call deflt(2,iv,liv,lv,v)
748 c$$$* 12 means fresh start, dont call deflt
750 c$$$* max num of fun calls
751 c$$$ if (maxfun.eq.0) maxfun=500
753 c$$$* max num of iterations
754 c$$$ if (maxmin.eq.0) maxmin=1000
756 c$$$* controls output
758 c$$$* selects output unit
760 c$$$c iv(21)=iout ! DEBUG
761 c$$$c iv(21)=8 ! DEBUG
762 c$$$* 1 means to print out result
764 c$$$c iv(22)=1 ! DEBUG
765 c$$$* 1 means to print out summary stats
767 c$$$c iv(23)=1 ! DEBUG
768 c$$$* 1 means to print initial x and d
770 c$$$c iv(24)=1 ! DEBUG
771 c$$$* min val for v(radfac) default is 0.1
773 c$$$* max val for v(radfac) default is 4.0
776 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
777 c$$$* the sumsl default is 0.1
779 c$$$* false conv if (act fnctn decrease) .lt. v(34)
780 c$$$* the sumsl default is 100*machep
781 c$$$ v(34)=v(34)/100.0D0
782 c$$$* absolute convergence
783 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
785 c$$$* relative convergence
786 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
788 c$$$* controls initial step size
790 c$$$* large vals of d correspond to small components of step
794 c$$$ do i=nphi+1,nvar
798 c$$$ call geom_to_var(nvar,x)
799 c$$$ IF (mask_r) THEN
800 c$$$ do i=1,nres ! Just in case...
804 c$$$ call x2xx(x,xx,nvar_restr)
805 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
806 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
809 c$$$c When minimizing ALL side-chains, etotal_sc is a little
810 c$$$c faster if we don't set mask_r
816 c$$$ call x2xx(x,xx,nvar_restr)
817 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
818 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
821 c$$$ call var_to_geom(nvar,x)
822 c$$$ call chainbuild_sc
829 c$$$C--------------------------------------------------------------------------
831 c$$$ subroutine chainbuild_sc
833 c$$$ include 'DIMENSIONS'
834 c$$$ include 'COMMON.VAR'
835 c$$$ include 'COMMON.INTERACT'
837 c$$$c Local variables
842 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
843 c$$$ call locate_side_chain(i)
850 c$$$C--------------------------------------------------------------------------
852 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
856 c$$$ include 'DIMENSIONS'
857 c$$$ include 'COMMON.DERIV'
858 c$$$ include 'COMMON.VAR'
859 c$$$ include 'COMMON.MINIM'
860 c$$$ include 'COMMON.IOUNITS'
862 c$$$c Input arguments
864 c$$$ double precision x(maxvar)
865 c$$$ double precision ufparm
868 c$$$c Input/Output arguments
870 c$$$ integer uiparm(1)
871 c$$$ double precision urparm(1)
873 c$$$c Output arguments
874 c$$$ double precision f
876 c$$$c Local variables
877 c$$$ double precision energia(0:n_ene)
879 c$$$c Variables used to intercept NaNs
880 c$$$ double precision x_sum
889 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
892 c$$$ x_sum=x_sum+x(i_NAN)
894 c$$$c Calculate the energy only if the coordinates are ok
895 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
896 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
902 c$$$ call var_to_geom_restr(n,x)
904 c$$$ call chainbuild_sc
905 c$$$ call etotal_sc(energia(0))
907 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
916 c$$$c-------------------------------------------------------
918 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
922 c$$$ include 'DIMENSIONS'
923 c$$$ include 'COMMON.CHAIN'
924 c$$$ include 'COMMON.DERIV'
925 c$$$ include 'COMMON.VAR'
926 c$$$ include 'COMMON.INTERACT'
927 c$$$ include 'COMMON.MINIM'
929 c$$$c Input arguments
931 c$$$ double precision x(maxvar)
932 c$$$ double precision ufparm
935 c$$$c Input/Output arguments
937 c$$$ integer uiparm(1)
938 c$$$ double precision urparm(1)
940 c$$$c Output arguments
941 c$$$ double precision g(maxvar)
943 c$$$c Local variables
944 c$$$ double precision f,gphii,gthetai,galphai,gomegai
945 c$$$ integer ig,ind,i,j,k,igall,ij
949 c$$$ if (nf-nfl+1) 20,30,40
950 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
951 c$$$c write (iout,*) 'grad 20'
952 c$$$ if (nf.eq.0) return
954 c$$$ 30 call var_to_geom_restr(n,x)
955 c$$$ call chainbuild_sc
957 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
961 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
967 c$$$ IF (mask_phi(i+2).eq.1) THEN
972 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
973 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
979 c$$$ ind=ind+nres-1-i
986 c$$$ IF (mask_theta(i+2).eq.1) THEN
992 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
993 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
998 c$$$ ind=ind+nres-1-i
1003 c$$$ if (itype(i).ne.10) then
1004 c$$$ IF (mask_side(i).eq.1) THEN
1008 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1017 c$$$ if (itype(i).ne.10) then
1018 c$$$ IF (mask_side(i).eq.1) THEN
1022 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1030 c$$$C Add the components corresponding to local energy terms.
1037 c$$$ if (mask_phi(i).eq.1) then
1039 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1045 c$$$ if (mask_theta(i).eq.1) then
1047 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1053 c$$$ if (itype(i).ne.10) then
1055 c$$$ if (mask_side(i).eq.1) then
1057 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1064 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1070 c$$$C-----------------------------------------------------------------------------
1072 c$$$ subroutine etotal_sc(energy_sc)
1076 c$$$ include 'DIMENSIONS'
1077 c$$$ include 'COMMON.VAR'
1078 c$$$ include 'COMMON.INTERACT'
1079 c$$$ include 'COMMON.DERIV'
1080 c$$$ include 'COMMON.FFIELD'
1082 c$$$c Output arguments
1083 c$$$ double precision energy_sc(0:n_ene)
1085 c$$$c Local variables
1086 c$$$ double precision evdw,escloc
1091 c$$$ energy_sc(i)=0.0D0
1094 c$$$ if (mask_r) then
1095 c$$$ call egb_sc(evdw)
1096 c$$$ call esc_sc(escloc)
1099 c$$$ call esc(escloc)
1102 c$$$ if (evdw.eq.1.0D20) then
1103 c$$$ energy_sc(0)=evdw
1105 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1107 c$$$ energy_sc(1)=evdw
1108 c$$$ energy_sc(12)=escloc
1111 c$$$C Sum up the components of the Cartesian gradient.
1115 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1122 c$$$C-----------------------------------------------------------------------------
1124 c$$$ subroutine egb_sc(evdw)
1126 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1127 c$$$C assuming the Gay-Berne potential of interaction.
1129 c$$$ implicit real*8 (a-h,o-z)
1130 c$$$ include 'DIMENSIONS'
1131 c$$$ include 'COMMON.GEO'
1132 c$$$ include 'COMMON.VAR'
1133 c$$$ include 'COMMON.LOCAL'
1134 c$$$ include 'COMMON.CHAIN'
1135 c$$$ include 'COMMON.DERIV'
1136 c$$$ include 'COMMON.NAMES'
1137 c$$$ include 'COMMON.INTERACT'
1138 c$$$ include 'COMMON.IOUNITS'
1139 c$$$ include 'COMMON.CALC'
1140 c$$$ include 'COMMON.CONTROL'
1143 c$$$ energy_dec=.false.
1144 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1147 c$$$c if (icall.eq.0) lprn=.false.
1149 c$$$ do i=iatsc_s,iatsc_e
1151 c$$$ itypi1=itype(i+1)
1155 c$$$ dxi=dc_norm(1,nres+i)
1156 c$$$ dyi=dc_norm(2,nres+i)
1157 c$$$ dzi=dc_norm(3,nres+i)
1158 c$$$c dsci_inv=dsc_inv(itypi)
1159 c$$$ dsci_inv=vbld_inv(i+nres)
1160 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1161 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1163 c$$$C Calculate SC interaction energy.
1165 c$$$ do iint=1,nint_gr(i)
1166 c$$$ do j=istart(i,iint),iend(i,iint)
1167 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1170 c$$$c dscj_inv=dsc_inv(itypj)
1171 c$$$ dscj_inv=vbld_inv(j+nres)
1172 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1173 c$$$c & 1.0d0/vbld(j+nres)
1174 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1175 c$$$ sig0ij=sigma(itypi,itypj)
1176 c$$$ chi1=chi(itypi,itypj)
1177 c$$$ chi2=chi(itypj,itypi)
1178 c$$$ chi12=chi1*chi2
1179 c$$$ chip1=chip(itypi)
1180 c$$$ chip2=chip(itypj)
1181 c$$$ chip12=chip1*chip2
1182 c$$$ alf1=alp(itypi)
1183 c$$$ alf2=alp(itypj)
1184 c$$$ alf12=0.5D0*(alf1+alf2)
1185 c$$$C For diagnostics only!!!
1195 c$$$ xj=c(1,nres+j)-xi
1196 c$$$ yj=c(2,nres+j)-yi
1197 c$$$ zj=c(3,nres+j)-zi
1198 c$$$ dxj=dc_norm(1,nres+j)
1199 c$$$ dyj=dc_norm(2,nres+j)
1200 c$$$ dzj=dc_norm(3,nres+j)
1201 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1202 c$$$c write (iout,*) "j",j," dc_norm",
1203 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1204 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1205 c$$$ rij=dsqrt(rrij)
1206 c$$$C Calculate angle-dependent terms of energy and contributions to their
1208 c$$$ call sc_angular
1209 c$$$ sigsq=1.0D0/sigsq
1210 c$$$ sig=sig0ij*dsqrt(sigsq)
1211 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1212 c$$$c for diagnostics; uncomment
1213 c$$$c rij_shift=1.2*sig0ij
1214 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1215 c$$$ if (rij_shift.le.0.0D0) then
1217 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1218 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1219 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1222 c$$$ sigder=-sig*sigsq
1223 c$$$c---------------------------------------------------------------
1224 c$$$ rij_shift=1.0D0/rij_shift
1225 c$$$ fac=rij_shift**expon
1226 c$$$ e1=fac*fac*aa(itypi,itypj)
1227 c$$$ e2=fac*bb(itypi,itypj)
1228 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1229 c$$$ eps2der=evdwij*eps3rt
1230 c$$$ eps3der=evdwij*eps2rt
1231 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1232 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1233 c$$$ evdwij=evdwij*eps2rt*eps3rt
1234 c$$$ evdw=evdw+evdwij
1236 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1237 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1238 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1239 c$$$ & restyp(itypi),i,restyp(itypj),j,
1240 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1241 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1242 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1246 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1247 c$$$ & 'evdw',i,j,evdwij
1249 c$$$C Calculate gradient components.
1250 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1251 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1252 c$$$ sigder=fac*sigder
1255 c$$$C Calculate the radial part of the gradient
1259 c$$$C Calculate angular part of the gradient.
1265 c$$$ energy_dec=.false.
1269 c$$$c-----------------------------------------------------------------------------
1271 c$$$ subroutine esc_sc(escloc)
1272 c$$$C Calculate the local energy of a side chain and its derivatives in the
1273 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1274 c$$$C ALPHA and OMEGA.
1275 c$$$ implicit real*8 (a-h,o-z)
1276 c$$$ include 'DIMENSIONS'
1277 c$$$ include 'COMMON.GEO'
1278 c$$$ include 'COMMON.LOCAL'
1279 c$$$ include 'COMMON.VAR'
1280 c$$$ include 'COMMON.INTERACT'
1281 c$$$ include 'COMMON.DERIV'
1282 c$$$ include 'COMMON.CHAIN'
1283 c$$$ include 'COMMON.IOUNITS'
1284 c$$$ include 'COMMON.NAMES'
1285 c$$$ include 'COMMON.FFIELD'
1286 c$$$ include 'COMMON.CONTROL'
1287 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1288 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1289 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1290 c$$$ delta=0.02d0*pi
1292 c$$$c write (iout,'(a)') 'ESC'
1293 c$$$ do i=loc_start,loc_end
1294 c$$$ IF (mask_side(i).eq.1) THEN
1296 c$$$ if (it.eq.10) goto 1
1297 c$$$ nlobit=nlob(it)
1298 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1299 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1300 c$$$ theti=theta(i+1)-pipol
1301 c$$$ x(1)=dtan(theti)
1305 c$$$ if (x(2).gt.pi-delta) then
1307 c$$$ xtemp(2)=pi-delta
1309 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1311 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1312 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1313 c$$$ & escloci,dersc(2))
1314 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1315 c$$$ & ddersc0(1),dersc(1))
1316 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1317 c$$$ & ddersc0(3),dersc(3))
1318 c$$$ xtemp(2)=pi-delta
1319 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1321 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1322 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1323 c$$$ & dersc0(2),esclocbi,dersc02)
1324 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1325 c$$$ & dersc12,dersc01)
1326 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1327 c$$$ dersc0(1)=dersc01
1328 c$$$ dersc0(2)=dersc02
1329 c$$$ dersc0(3)=0.0d0
1331 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1333 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1334 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1335 c$$$c & esclocbi,ss,ssd
1336 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1337 c$$$c escloci=esclocbi
1338 c$$$c write (iout,*) escloci
1339 c$$$ else if (x(2).lt.delta) then
1343 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1345 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1346 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1347 c$$$ & escloci,dersc(2))
1348 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1349 c$$$ & ddersc0(1),dersc(1))
1350 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1351 c$$$ & ddersc0(3),dersc(3))
1353 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1355 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1356 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1357 c$$$ & dersc0(2),esclocbi,dersc02)
1358 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1359 c$$$ & dersc12,dersc01)
1360 c$$$ dersc0(1)=dersc01
1361 c$$$ dersc0(2)=dersc02
1362 c$$$ dersc0(3)=0.0d0
1363 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1365 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1367 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1368 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1369 c$$$c & esclocbi,ss,ssd
1370 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1371 c$$$c write (iout,*) escloci
1373 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1376 c$$$ escloc=escloc+escloci
1377 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1378 c$$$ & 'escloc',i,escloci
1379 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1381 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1382 c$$$ & wscloc*dersc(1)
1383 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1384 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1391 c$$$C-----------------------------------------------------------------------------
1393 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1395 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1396 c$$$C assuming the Gay-Berne potential of interaction.
1398 c$$$ implicit real*8 (a-h,o-z)
1399 c$$$ include 'DIMENSIONS'
1400 c$$$ include 'COMMON.GEO'
1401 c$$$ include 'COMMON.VAR'
1402 c$$$ include 'COMMON.LOCAL'
1403 c$$$ include 'COMMON.CHAIN'
1404 c$$$ include 'COMMON.DERIV'
1405 c$$$ include 'COMMON.NAMES'
1406 c$$$ include 'COMMON.INTERACT'
1407 c$$$ include 'COMMON.IOUNITS'
1408 c$$$ include 'COMMON.CALC'
1409 c$$$ include 'COMMON.CONTROL'
1412 c$$$ energy_dec=.false.
1413 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1417 c$$$c$$$ do i=iatsc_s,iatsc_e
1420 c$$$ itypi1=itype(i+1)
1424 c$$$ dxi=dc_norm(1,nres+i)
1425 c$$$ dyi=dc_norm(2,nres+i)
1426 c$$$ dzi=dc_norm(3,nres+i)
1427 c$$$c dsci_inv=dsc_inv(itypi)
1428 c$$$ dsci_inv=vbld_inv(i+nres)
1429 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1430 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1432 c$$$C Calculate SC interaction energy.
1434 c$$$c$$$ do iint=1,nint_gr(i)
1435 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1439 c$$$c dscj_inv=dsc_inv(itypj)
1440 c$$$ dscj_inv=vbld_inv(j+nres)
1441 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1442 c$$$c & 1.0d0/vbld(j+nres)
1443 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1444 c$$$ sig0ij=sigma(itypi,itypj)
1445 c$$$ chi1=chi(itypi,itypj)
1446 c$$$ chi2=chi(itypj,itypi)
1447 c$$$ chi12=chi1*chi2
1448 c$$$ chip1=chip(itypi)
1449 c$$$ chip2=chip(itypj)
1450 c$$$ chip12=chip1*chip2
1451 c$$$ alf1=alp(itypi)
1452 c$$$ alf2=alp(itypj)
1453 c$$$ alf12=0.5D0*(alf1+alf2)
1454 c$$$C For diagnostics only!!!
1464 c$$$ xj=c(1,nres+j)-xi
1465 c$$$ yj=c(2,nres+j)-yi
1466 c$$$ zj=c(3,nres+j)-zi
1467 c$$$ dxj=dc_norm(1,nres+j)
1468 c$$$ dyj=dc_norm(2,nres+j)
1469 c$$$ dzj=dc_norm(3,nres+j)
1470 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1471 c$$$c write (iout,*) "j",j," dc_norm",
1472 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1473 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1474 c$$$ rij=dsqrt(rrij)
1475 c$$$C Calculate angle-dependent terms of energy and contributions to their
1477 c$$$ call sc_angular
1478 c$$$ sigsq=1.0D0/sigsq
1479 c$$$ sig=sig0ij*dsqrt(sigsq)
1480 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1481 c$$$c for diagnostics; uncomment
1482 c$$$c rij_shift=1.2*sig0ij
1483 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1484 c$$$ if (rij_shift.le.0.0D0) then
1486 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1487 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1488 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1491 c$$$ sigder=-sig*sigsq
1492 c$$$c---------------------------------------------------------------
1493 c$$$ rij_shift=1.0D0/rij_shift
1494 c$$$ fac=rij_shift**expon
1495 c$$$ e1=fac*fac*aa(itypi,itypj)
1496 c$$$ e2=fac*bb(itypi,itypj)
1497 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1498 c$$$ eps2der=evdwij*eps3rt
1499 c$$$ eps3der=evdwij*eps2rt
1500 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1501 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1502 c$$$ evdwij=evdwij*eps2rt*eps3rt
1503 c$$$ evdw=evdw+evdwij
1505 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1506 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1507 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1508 c$$$ & restyp(itypi),i,restyp(itypj),j,
1509 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1510 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1511 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1515 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1516 c$$$ & 'evdw',i,j,evdwij
1518 c$$$C Calculate gradient components.
1519 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1520 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1521 c$$$ sigder=fac*sigder
1524 c$$$C Calculate the radial part of the gradient
1528 c$$$C Calculate angular part of the gradient.
1531 c$$$c$$$ enddo ! iint
1533 c$$$ energy_dec=.false.
1537 c$$$C-----------------------------------------------------------------------------
1539 c$$$ subroutine perturb_side_chain(i,angle)
1543 c$$$ include 'DIMENSIONS'
1544 c$$$ include 'COMMON.CHAIN'
1545 c$$$ include 'COMMON.GEO'
1546 c$$$ include 'COMMON.VAR'
1547 c$$$ include 'COMMON.LOCAL'
1548 c$$$ include 'COMMON.IOUNITS'
1550 c$$$c External functions
1551 c$$$ external ran_number
1552 c$$$ double precision ran_number
1554 c$$$c Input arguments
1556 c$$$ double precision angle ! In degrees
1558 c$$$c Local variables
1560 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1564 c$$$ rad_ang=angle*deg2rad
1567 c$$$ do while (length.lt.0.01)
1568 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1569 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1570 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1571 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1572 c$$$ + rand_v(3)*rand_v(3)
1573 c$$$ length=sqrt(length)
1574 c$$$ rand_v(1)=rand_v(1)/length
1575 c$$$ rand_v(2)=rand_v(2)/length
1576 c$$$ rand_v(3)=rand_v(3)/length
1577 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1578 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1579 c$$$ length=1.0D0-cost*cost
1580 c$$$ if (length.lt.0.0D0) length=0.0D0
1581 c$$$ length=sqrt(length)
1582 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1583 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1584 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1586 c$$$ rand_v(1)=rand_v(1)/length
1587 c$$$ rand_v(2)=rand_v(2)/length
1588 c$$$ rand_v(3)=rand_v(3)/length
1590 c$$$ cost=dcos(rad_ang)
1591 c$$$ sint=dsin(rad_ang)
1592 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1593 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1594 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1595 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1596 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1597 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1598 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1599 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1600 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1602 c$$$ call chainbuild_cart
1607 c$$$c----------------------------------------------------------------------------
1609 c$$$ subroutine ss_relax3(i_in,j_in)
1613 c$$$ include 'DIMENSIONS'
1614 c$$$ include 'COMMON.VAR'
1615 c$$$ include 'COMMON.CHAIN'
1616 c$$$ include 'COMMON.IOUNITS'
1617 c$$$ include 'COMMON.INTERACT'
1619 c$$$c External functions
1620 c$$$ external ran_number
1621 c$$$ double precision ran_number
1623 c$$$c Input arguments
1624 c$$$ integer i_in,j_in
1626 c$$$c Local variables
1627 c$$$ double precision energy_sc(0:n_ene),etot
1628 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1629 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1630 c$$$ integer n,i_pert,i
1631 c$$$ logical notdone
1640 c$$$ mask_side(i_in)=1
1641 c$$$ mask_side(j_in)=1
1643 c$$$ call etotal_sc(energy_sc)
1644 c$$$ etot=energy_sc(0)
1645 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1646 c$$$c + energy_sc(1),energy_sc(12)
1650 c$$$ do while (notdone)
1651 c$$$ if (mod(n,2).eq.0) then
1659 c$$$ org_dc(i)=dc(i,i_pert+nres)
1660 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1661 c$$$ org_c(i)=c(i,i_pert+nres)
1663 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1664 c$$$ call perturb_side_chain(i_pert,ang_pert)
1665 c$$$ call etotal_sc(energy_sc)
1666 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1667 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1668 c$$$ if (rand_fact.lt.exp_fact) then
1669 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1670 c$$$c + energy_sc(1),energy_sc(12)
1671 c$$$ etot=energy_sc(0)
1673 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1674 c$$$c + energy_sc(1),energy_sc(12)
1676 c$$$ dc(i,i_pert+nres)=org_dc(i)
1677 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1678 c$$$ c(i,i_pert+nres)=org_c(i)
1682 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1690 c$$$c----------------------------------------------------------------------------
1692 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1694 c$$$ include 'DIMENSIONS'
1696 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1697 c$$$*********************************************************************
1698 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1699 c$$$* the calling subprogram. *
1700 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1701 c$$$* calculated in the usual pythagorean way. *
1702 c$$$* absolute convergence occurs when the function is within v(31) of *
1703 c$$$* zero. unless you know the minimum value in advance, abs convg *
1704 c$$$* is probably not useful. *
1705 c$$$* relative convergence is when the model predicts that the function *
1706 c$$$* will decrease by less than v(32)*abs(fun). *
1707 c$$$*********************************************************************
1708 c$$$ include 'COMMON.IOUNITS'
1709 c$$$ include 'COMMON.VAR'
1710 c$$$ include 'COMMON.GEO'
1711 c$$$ include 'COMMON.MINIM'
1712 c$$$ include 'COMMON.CHAIN'
1714 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1715 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1716 c$$$ + orig_ss_dist(maxres2,maxres2)
1718 c$$$ double precision etot
1719 c$$$ integer iretcode,nfun,i_in,j_in
1722 c$$$ double precision dist
1723 c$$$ external ss_func,fdum
1724 c$$$ double precision ss_func,fdum
1726 c$$$ integer iv(liv),uiparm(2)
1727 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1731 c$$$ call deflt(2,iv,liv,lv,v)
1732 c$$$* 12 means fresh start, dont call deflt
1734 c$$$* max num of fun calls
1735 c$$$ if (maxfun.eq.0) maxfun=500
1737 c$$$* max num of iterations
1738 c$$$ if (maxmin.eq.0) maxmin=1000
1740 c$$$* controls output
1742 c$$$* selects output unit
1745 c$$$* 1 means to print out result
1747 c$$$* 1 means to print out summary stats
1749 c$$$* 1 means to print initial x and d
1751 c$$$* min val for v(radfac) default is 0.1
1753 c$$$* max val for v(radfac) default is 4.0
1756 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1757 c$$$* the sumsl default is 0.1
1759 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1760 c$$$* the sumsl default is 100*machep
1761 c$$$ v(34)=v(34)/100.0D0
1762 c$$$* absolute convergence
1763 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1766 c$$$* relative convergence
1767 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1770 c$$$* controls initial step size
1772 c$$$* large vals of d correspond to small components of step
1779 c$$$ orig_ss_dc(j,i)=dc(j,i)
1782 c$$$ call geom_to_var(nvar,orig_ss_var)
1786 c$$$ orig_ss_dist(j,i)=dist(j,i)
1787 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1788 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1789 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1801 c$$$ if (ialph(i,1).gt.0) then
1804 c$$$ x(k)=dc(j,i+nres)
1811 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1814 c$$$ nfun=iv(6)+iv(30)
1824 c$$$ if (ialph(i,1).gt.0) then
1827 c$$$ dc(j,i+nres)=x(k)
1831 c$$$ call chainbuild_cart
1836 c$$$C-----------------------------------------------------------------------------
1838 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1840 c$$$ include 'DIMENSIONS'
1841 c$$$ include 'COMMON.DERIV'
1842 c$$$ include 'COMMON.IOUNITS'
1843 c$$$ include 'COMMON.VAR'
1844 c$$$ include 'COMMON.CHAIN'
1845 c$$$ include 'COMMON.INTERACT'
1846 c$$$ include 'COMMON.SBRIDGE'
1848 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1849 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1850 c$$$ + orig_ss_dist(maxres2,maxres2)
1853 c$$$ double precision x(maxres6)
1855 c$$$ double precision f
1856 c$$$ integer uiparm(2)
1857 c$$$ real*8 urparm(1)
1858 c$$$ external ufparm
1859 c$$$ double precision ufparm
1862 c$$$ double precision dist
1864 c$$$ integer i,j,k,ss_i,ss_j
1865 c$$$ double precision tempf,var(maxvar)
1880 c$$$ if (ialph(i,1).gt.0) then
1883 c$$$ dc(j,i+nres)=x(k)
1887 c$$$ call chainbuild_cart
1889 c$$$ call geom_to_var(nvar,var)
1891 c$$$c Constraints on all angles
1893 c$$$ tempf=var(i)-orig_ss_var(i)
1894 c$$$ f=f+tempf*tempf
1897 c$$$c Constraints on all distances
1899 c$$$ if (i.gt.1) then
1900 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1901 c$$$ f=f+tempf*tempf
1904 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1905 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1906 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1907 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1908 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1909 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1910 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1911 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1915 c$$$c Constraints for the relevant CYS-CYS
1916 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1917 c$$$ f=f+tempf*tempf
1918 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
1920 c$$$c$$$ if (nf.ne.nfl) then
1921 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1922 c$$$c$$$ + f,dist(5+nres,14+nres)
1930 c$$$C-----------------------------------------------------------------------------