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
518 include 'COMMON.SBRIDGE'
519 include 'COMMON.CHAIN'
520 include 'COMMON.IOUNITS'
521 include 'COMMON.SETUP'
529 double precision emin
531 integer diff,allflag(maxdim),allnss,
532 & allihpb(maxdim),alljhpb(maxdim),
533 & newnss,newihpb(maxdim),newjhpb(maxdim)
535 integer i_newnss(max_fg_procs),displ(max_fg_procs)
536 integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
541 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
550 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
554 if (allflag(i).eq.0 .and.
555 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
556 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
560 if (emin.lt.1.0d300) then
563 if (allflag(i).eq.0 .and.
564 & (allihpb(i).eq.allihpb(imin) .or.
565 & alljhpb(i).eq.allihpb(imin) .or.
566 & allihpb(i).eq.alljhpb(imin) .or.
567 & alljhpb(i).eq.alljhpb(imin))) then
574 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
578 if (allflag(i).eq.1) then
580 newihpb(newnss)=allihpb(i)
581 newjhpb(newnss)=alljhpb(i)
586 if (nfgtasks.gt.1)then
588 call MPI_Reduce(newnss,g_newnss,1,
589 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
590 call MPI_Gather(newnss,1,MPI_INTEGER,
591 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
594 displ(i)=i_newnss(i-1)+displ(i-1)
596 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
597 & g_newihpb,i_newnss,displ,MPI_INTEGER,
599 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
600 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
602 if(fg_rank.eq.0) then
603 c print *,'g_newnss',g_newnss
604 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
605 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
608 newihpb(i)=g_newihpb(i)
609 newjhpb(i)=g_newjhpb(i)
617 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
622 if (idssb(i).eq.newihpb(j) .and.
623 & jdssb(i).eq.newjhpb(j)) found=.true.
627 if (.not.found.and.fg_rank.eq.0)
628 & write(iout,'(a15,f12.2,f8.1,2i5)')
629 & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
637 if (newihpb(i).eq.idssb(j) .and.
638 & newjhpb(i).eq.jdssb(j)) found=.true.
642 if (.not.found.and.fg_rank.eq.0)
643 & write(iout,'(a15,f12.2,f8.1,2i5)')
644 & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
659 c$$$c-----------------------------------------------------------------------------
661 c$$$ subroutine ss_relax(i_in,j_in)
665 c$$$ include 'DIMENSIONS'
666 c$$$ include 'COMMON.VAR'
667 c$$$ include 'COMMON.CHAIN'
668 c$$$ include 'COMMON.IOUNITS'
669 c$$$ include 'COMMON.INTERACT'
671 c$$$c Input arguments
672 c$$$ integer i_in,j_in
674 c$$$c Local variables
675 c$$$ integer i,iretcode,nfun_sc
677 c$$$ double precision var(maxvar),e_sc,etot
684 c$$$ mask_side(i_in)=1
685 c$$$ mask_side(j_in)=1
687 c$$$c Minimize the two selected side-chains
688 c$$$ call overlap_sc(scfail) ! Better not fail!
689 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
696 c$$$c-------------------------------------------------------------
698 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
699 c$$$c Minimize side-chains only, starting from geom but without modifying
701 c$$$c If mask_r is already set, only the selected side-chains are minimized,
702 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
706 c$$$ include 'DIMENSIONS'
707 c$$$ include 'COMMON.IOUNITS'
708 c$$$ include 'COMMON.VAR'
709 c$$$ include 'COMMON.CHAIN'
710 c$$$ include 'COMMON.GEO'
711 c$$$ include 'COMMON.MINIM'
713 c$$$ common /srutu/ icall
715 c$$$c Output arguments
716 c$$$ double precision etot_sc
717 c$$$ integer iretcode,nfun
719 c$$$c External functions/subroutines
720 c$$$ external func_sc,grad_sc,fdum
722 c$$$c Local variables
724 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
726 c$$$ double precision rdum(1)
727 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
729 c$$$ integer i,nvar_restr
732 c$$$cmc start_minim=.true.
733 c$$$ call deflt(2,iv,liv,lv,v)
734 c$$$* 12 means fresh start, dont call deflt
736 c$$$* max num of fun calls
737 c$$$ if (maxfun.eq.0) maxfun=500
739 c$$$* max num of iterations
740 c$$$ if (maxmin.eq.0) maxmin=1000
742 c$$$* controls output
744 c$$$* selects output unit
746 c$$$c iv(21)=iout ! DEBUG
747 c$$$c iv(21)=8 ! DEBUG
748 c$$$* 1 means to print out result
750 c$$$c iv(22)=1 ! DEBUG
751 c$$$* 1 means to print out summary stats
753 c$$$c iv(23)=1 ! DEBUG
754 c$$$* 1 means to print initial x and d
756 c$$$c iv(24)=1 ! DEBUG
757 c$$$* min val for v(radfac) default is 0.1
759 c$$$* max val for v(radfac) default is 4.0
762 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
763 c$$$* the sumsl default is 0.1
765 c$$$* false conv if (act fnctn decrease) .lt. v(34)
766 c$$$* the sumsl default is 100*machep
767 c$$$ v(34)=v(34)/100.0D0
768 c$$$* absolute convergence
769 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
771 c$$$* relative convergence
772 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
774 c$$$* controls initial step size
776 c$$$* large vals of d correspond to small components of step
780 c$$$ do i=nphi+1,nvar
784 c$$$ call geom_to_var(nvar,x)
785 c$$$ IF (mask_r) THEN
786 c$$$ do i=1,nres ! Just in case...
790 c$$$ call x2xx(x,xx,nvar_restr)
791 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
792 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
795 c$$$c When minimizing ALL side-chains, etotal_sc is a little
796 c$$$c faster if we don't set mask_r
802 c$$$ call x2xx(x,xx,nvar_restr)
803 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
804 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
807 c$$$ call var_to_geom(nvar,x)
808 c$$$ call chainbuild_sc
815 c$$$C--------------------------------------------------------------------------
817 c$$$ subroutine chainbuild_sc
819 c$$$ include 'DIMENSIONS'
820 c$$$ include 'COMMON.VAR'
821 c$$$ include 'COMMON.INTERACT'
823 c$$$c Local variables
828 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
829 c$$$ call locate_side_chain(i)
836 c$$$C--------------------------------------------------------------------------
838 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
842 c$$$ include 'DIMENSIONS'
843 c$$$ include 'COMMON.DERIV'
844 c$$$ include 'COMMON.VAR'
845 c$$$ include 'COMMON.MINIM'
846 c$$$ include 'COMMON.IOUNITS'
848 c$$$c Input arguments
850 c$$$ double precision x(maxvar)
851 c$$$ double precision ufparm
854 c$$$c Input/Output arguments
856 c$$$ integer uiparm(1)
857 c$$$ double precision urparm(1)
859 c$$$c Output arguments
860 c$$$ double precision f
862 c$$$c Local variables
863 c$$$ double precision energia(0:n_ene)
865 c$$$c Variables used to intercept NaNs
866 c$$$ double precision x_sum
875 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
878 c$$$ x_sum=x_sum+x(i_NAN)
880 c$$$c Calculate the energy only if the coordinates are ok
881 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
882 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
888 c$$$ call var_to_geom_restr(n,x)
890 c$$$ call chainbuild_sc
891 c$$$ call etotal_sc(energia(0))
893 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
902 c$$$c-------------------------------------------------------
904 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
908 c$$$ include 'DIMENSIONS'
909 c$$$ include 'COMMON.CHAIN'
910 c$$$ include 'COMMON.DERIV'
911 c$$$ include 'COMMON.VAR'
912 c$$$ include 'COMMON.INTERACT'
913 c$$$ include 'COMMON.MINIM'
915 c$$$c Input arguments
917 c$$$ double precision x(maxvar)
918 c$$$ double precision ufparm
921 c$$$c Input/Output arguments
923 c$$$ integer uiparm(1)
924 c$$$ double precision urparm(1)
926 c$$$c Output arguments
927 c$$$ double precision g(maxvar)
929 c$$$c Local variables
930 c$$$ double precision f,gphii,gthetai,galphai,gomegai
931 c$$$ integer ig,ind,i,j,k,igall,ij
935 c$$$ if (nf-nfl+1) 20,30,40
936 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
937 c$$$c write (iout,*) 'grad 20'
938 c$$$ if (nf.eq.0) return
940 c$$$ 30 call var_to_geom_restr(n,x)
941 c$$$ call chainbuild_sc
943 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
947 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
953 c$$$ IF (mask_phi(i+2).eq.1) THEN
958 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
959 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
965 c$$$ ind=ind+nres-1-i
972 c$$$ IF (mask_theta(i+2).eq.1) THEN
978 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
979 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
984 c$$$ ind=ind+nres-1-i
989 c$$$ if (itype(i).ne.10) then
990 c$$$ IF (mask_side(i).eq.1) THEN
994 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1003 c$$$ if (itype(i).ne.10) then
1004 c$$$ IF (mask_side(i).eq.1) THEN
1008 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1016 c$$$C Add the components corresponding to local energy terms.
1023 c$$$ if (mask_phi(i).eq.1) then
1025 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1031 c$$$ if (mask_theta(i).eq.1) then
1033 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1039 c$$$ if (itype(i).ne.10) then
1041 c$$$ if (mask_side(i).eq.1) then
1043 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1050 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1056 c$$$C-----------------------------------------------------------------------------
1058 c$$$ subroutine etotal_sc(energy_sc)
1062 c$$$ include 'DIMENSIONS'
1063 c$$$ include 'COMMON.VAR'
1064 c$$$ include 'COMMON.INTERACT'
1065 c$$$ include 'COMMON.DERIV'
1066 c$$$ include 'COMMON.FFIELD'
1068 c$$$c Output arguments
1069 c$$$ double precision energy_sc(0:n_ene)
1071 c$$$c Local variables
1072 c$$$ double precision evdw,escloc
1077 c$$$ energy_sc(i)=0.0D0
1080 c$$$ if (mask_r) then
1081 c$$$ call egb_sc(evdw)
1082 c$$$ call esc_sc(escloc)
1085 c$$$ call esc(escloc)
1088 c$$$ if (evdw.eq.1.0D20) then
1089 c$$$ energy_sc(0)=evdw
1091 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1093 c$$$ energy_sc(1)=evdw
1094 c$$$ energy_sc(12)=escloc
1097 c$$$C Sum up the components of the Cartesian gradient.
1101 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1108 c$$$C-----------------------------------------------------------------------------
1110 c$$$ subroutine egb_sc(evdw)
1112 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1113 c$$$C assuming the Gay-Berne potential of interaction.
1115 c$$$ implicit real*8 (a-h,o-z)
1116 c$$$ include 'DIMENSIONS'
1117 c$$$ include 'COMMON.GEO'
1118 c$$$ include 'COMMON.VAR'
1119 c$$$ include 'COMMON.LOCAL'
1120 c$$$ include 'COMMON.CHAIN'
1121 c$$$ include 'COMMON.DERIV'
1122 c$$$ include 'COMMON.NAMES'
1123 c$$$ include 'COMMON.INTERACT'
1124 c$$$ include 'COMMON.IOUNITS'
1125 c$$$ include 'COMMON.CALC'
1126 c$$$ include 'COMMON.CONTROL'
1129 c$$$ energy_dec=.false.
1130 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1133 c$$$c if (icall.eq.0) lprn=.false.
1135 c$$$ do i=iatsc_s,iatsc_e
1137 c$$$ itypi1=itype(i+1)
1141 c$$$ dxi=dc_norm(1,nres+i)
1142 c$$$ dyi=dc_norm(2,nres+i)
1143 c$$$ dzi=dc_norm(3,nres+i)
1144 c$$$c dsci_inv=dsc_inv(itypi)
1145 c$$$ dsci_inv=vbld_inv(i+nres)
1146 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1147 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1149 c$$$C Calculate SC interaction energy.
1151 c$$$ do iint=1,nint_gr(i)
1152 c$$$ do j=istart(i,iint),iend(i,iint)
1153 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1156 c$$$c dscj_inv=dsc_inv(itypj)
1157 c$$$ dscj_inv=vbld_inv(j+nres)
1158 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1159 c$$$c & 1.0d0/vbld(j+nres)
1160 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1161 c$$$ sig0ij=sigma(itypi,itypj)
1162 c$$$ chi1=chi(itypi,itypj)
1163 c$$$ chi2=chi(itypj,itypi)
1164 c$$$ chi12=chi1*chi2
1165 c$$$ chip1=chip(itypi)
1166 c$$$ chip2=chip(itypj)
1167 c$$$ chip12=chip1*chip2
1168 c$$$ alf1=alp(itypi)
1169 c$$$ alf2=alp(itypj)
1170 c$$$ alf12=0.5D0*(alf1+alf2)
1171 c$$$C For diagnostics only!!!
1181 c$$$ xj=c(1,nres+j)-xi
1182 c$$$ yj=c(2,nres+j)-yi
1183 c$$$ zj=c(3,nres+j)-zi
1184 c$$$ dxj=dc_norm(1,nres+j)
1185 c$$$ dyj=dc_norm(2,nres+j)
1186 c$$$ dzj=dc_norm(3,nres+j)
1187 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1188 c$$$c write (iout,*) "j",j," dc_norm",
1189 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1190 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1191 c$$$ rij=dsqrt(rrij)
1192 c$$$C Calculate angle-dependent terms of energy and contributions to their
1194 c$$$ call sc_angular
1195 c$$$ sigsq=1.0D0/sigsq
1196 c$$$ sig=sig0ij*dsqrt(sigsq)
1197 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1198 c$$$c for diagnostics; uncomment
1199 c$$$c rij_shift=1.2*sig0ij
1200 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1201 c$$$ if (rij_shift.le.0.0D0) then
1203 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1204 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1205 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1208 c$$$ sigder=-sig*sigsq
1209 c$$$c---------------------------------------------------------------
1210 c$$$ rij_shift=1.0D0/rij_shift
1211 c$$$ fac=rij_shift**expon
1212 c$$$ e1=fac*fac*aa(itypi,itypj)
1213 c$$$ e2=fac*bb(itypi,itypj)
1214 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1215 c$$$ eps2der=evdwij*eps3rt
1216 c$$$ eps3der=evdwij*eps2rt
1217 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1218 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1219 c$$$ evdwij=evdwij*eps2rt*eps3rt
1220 c$$$ evdw=evdw+evdwij
1222 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1223 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1224 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1225 c$$$ & restyp(itypi),i,restyp(itypj),j,
1226 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1227 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1228 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1232 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1233 c$$$ & 'evdw',i,j,evdwij
1235 c$$$C Calculate gradient components.
1236 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1237 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1238 c$$$ sigder=fac*sigder
1241 c$$$C Calculate the radial part of the gradient
1245 c$$$C Calculate angular part of the gradient.
1251 c$$$ energy_dec=.false.
1255 c$$$c-----------------------------------------------------------------------------
1257 c$$$ subroutine esc_sc(escloc)
1258 c$$$C Calculate the local energy of a side chain and its derivatives in the
1259 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1260 c$$$C ALPHA and OMEGA.
1261 c$$$ implicit real*8 (a-h,o-z)
1262 c$$$ include 'DIMENSIONS'
1263 c$$$ include 'COMMON.GEO'
1264 c$$$ include 'COMMON.LOCAL'
1265 c$$$ include 'COMMON.VAR'
1266 c$$$ include 'COMMON.INTERACT'
1267 c$$$ include 'COMMON.DERIV'
1268 c$$$ include 'COMMON.CHAIN'
1269 c$$$ include 'COMMON.IOUNITS'
1270 c$$$ include 'COMMON.NAMES'
1271 c$$$ include 'COMMON.FFIELD'
1272 c$$$ include 'COMMON.CONTROL'
1273 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1274 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1275 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1276 c$$$ delta=0.02d0*pi
1278 c$$$c write (iout,'(a)') 'ESC'
1279 c$$$ do i=loc_start,loc_end
1280 c$$$ IF (mask_side(i).eq.1) THEN
1282 c$$$ if (it.eq.10) goto 1
1283 c$$$ nlobit=nlob(it)
1284 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1285 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1286 c$$$ theti=theta(i+1)-pipol
1287 c$$$ x(1)=dtan(theti)
1291 c$$$ if (x(2).gt.pi-delta) then
1293 c$$$ xtemp(2)=pi-delta
1295 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1297 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1298 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1299 c$$$ & escloci,dersc(2))
1300 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1301 c$$$ & ddersc0(1),dersc(1))
1302 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1303 c$$$ & ddersc0(3),dersc(3))
1304 c$$$ xtemp(2)=pi-delta
1305 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1307 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1308 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1309 c$$$ & dersc0(2),esclocbi,dersc02)
1310 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1311 c$$$ & dersc12,dersc01)
1312 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1313 c$$$ dersc0(1)=dersc01
1314 c$$$ dersc0(2)=dersc02
1315 c$$$ dersc0(3)=0.0d0
1317 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1319 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1320 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1321 c$$$c & esclocbi,ss,ssd
1322 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1323 c$$$c escloci=esclocbi
1324 c$$$c write (iout,*) escloci
1325 c$$$ else if (x(2).lt.delta) then
1329 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1331 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1332 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1333 c$$$ & escloci,dersc(2))
1334 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1335 c$$$ & ddersc0(1),dersc(1))
1336 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1337 c$$$ & ddersc0(3),dersc(3))
1339 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1341 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1342 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1343 c$$$ & dersc0(2),esclocbi,dersc02)
1344 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1345 c$$$ & dersc12,dersc01)
1346 c$$$ dersc0(1)=dersc01
1347 c$$$ dersc0(2)=dersc02
1348 c$$$ dersc0(3)=0.0d0
1349 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1351 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1353 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1354 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1355 c$$$c & esclocbi,ss,ssd
1356 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1357 c$$$c write (iout,*) escloci
1359 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1362 c$$$ escloc=escloc+escloci
1363 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1364 c$$$ & 'escloc',i,escloci
1365 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1367 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1368 c$$$ & wscloc*dersc(1)
1369 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1370 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1377 c$$$C-----------------------------------------------------------------------------
1379 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1381 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1382 c$$$C assuming the Gay-Berne potential of interaction.
1384 c$$$ implicit real*8 (a-h,o-z)
1385 c$$$ include 'DIMENSIONS'
1386 c$$$ include 'COMMON.GEO'
1387 c$$$ include 'COMMON.VAR'
1388 c$$$ include 'COMMON.LOCAL'
1389 c$$$ include 'COMMON.CHAIN'
1390 c$$$ include 'COMMON.DERIV'
1391 c$$$ include 'COMMON.NAMES'
1392 c$$$ include 'COMMON.INTERACT'
1393 c$$$ include 'COMMON.IOUNITS'
1394 c$$$ include 'COMMON.CALC'
1395 c$$$ include 'COMMON.CONTROL'
1398 c$$$ energy_dec=.false.
1399 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1403 c$$$c$$$ do i=iatsc_s,iatsc_e
1406 c$$$ itypi1=itype(i+1)
1410 c$$$ dxi=dc_norm(1,nres+i)
1411 c$$$ dyi=dc_norm(2,nres+i)
1412 c$$$ dzi=dc_norm(3,nres+i)
1413 c$$$c dsci_inv=dsc_inv(itypi)
1414 c$$$ dsci_inv=vbld_inv(i+nres)
1415 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1416 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1418 c$$$C Calculate SC interaction energy.
1420 c$$$c$$$ do iint=1,nint_gr(i)
1421 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1425 c$$$c dscj_inv=dsc_inv(itypj)
1426 c$$$ dscj_inv=vbld_inv(j+nres)
1427 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1428 c$$$c & 1.0d0/vbld(j+nres)
1429 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1430 c$$$ sig0ij=sigma(itypi,itypj)
1431 c$$$ chi1=chi(itypi,itypj)
1432 c$$$ chi2=chi(itypj,itypi)
1433 c$$$ chi12=chi1*chi2
1434 c$$$ chip1=chip(itypi)
1435 c$$$ chip2=chip(itypj)
1436 c$$$ chip12=chip1*chip2
1437 c$$$ alf1=alp(itypi)
1438 c$$$ alf2=alp(itypj)
1439 c$$$ alf12=0.5D0*(alf1+alf2)
1440 c$$$C For diagnostics only!!!
1450 c$$$ xj=c(1,nres+j)-xi
1451 c$$$ yj=c(2,nres+j)-yi
1452 c$$$ zj=c(3,nres+j)-zi
1453 c$$$ dxj=dc_norm(1,nres+j)
1454 c$$$ dyj=dc_norm(2,nres+j)
1455 c$$$ dzj=dc_norm(3,nres+j)
1456 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1457 c$$$c write (iout,*) "j",j," dc_norm",
1458 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1459 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1460 c$$$ rij=dsqrt(rrij)
1461 c$$$C Calculate angle-dependent terms of energy and contributions to their
1463 c$$$ call sc_angular
1464 c$$$ sigsq=1.0D0/sigsq
1465 c$$$ sig=sig0ij*dsqrt(sigsq)
1466 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1467 c$$$c for diagnostics; uncomment
1468 c$$$c rij_shift=1.2*sig0ij
1469 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1470 c$$$ if (rij_shift.le.0.0D0) then
1472 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1473 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1474 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1477 c$$$ sigder=-sig*sigsq
1478 c$$$c---------------------------------------------------------------
1479 c$$$ rij_shift=1.0D0/rij_shift
1480 c$$$ fac=rij_shift**expon
1481 c$$$ e1=fac*fac*aa(itypi,itypj)
1482 c$$$ e2=fac*bb(itypi,itypj)
1483 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1484 c$$$ eps2der=evdwij*eps3rt
1485 c$$$ eps3der=evdwij*eps2rt
1486 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1487 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1488 c$$$ evdwij=evdwij*eps2rt*eps3rt
1489 c$$$ evdw=evdw+evdwij
1491 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1492 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1493 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1494 c$$$ & restyp(itypi),i,restyp(itypj),j,
1495 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1496 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1497 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1501 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1502 c$$$ & 'evdw',i,j,evdwij
1504 c$$$C Calculate gradient components.
1505 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1506 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1507 c$$$ sigder=fac*sigder
1510 c$$$C Calculate the radial part of the gradient
1514 c$$$C Calculate angular part of the gradient.
1517 c$$$c$$$ enddo ! iint
1519 c$$$ energy_dec=.false.
1523 c$$$C-----------------------------------------------------------------------------
1525 c$$$ subroutine perturb_side_chain(i,angle)
1529 c$$$ include 'DIMENSIONS'
1530 c$$$ include 'COMMON.CHAIN'
1531 c$$$ include 'COMMON.GEO'
1532 c$$$ include 'COMMON.VAR'
1533 c$$$ include 'COMMON.LOCAL'
1534 c$$$ include 'COMMON.IOUNITS'
1536 c$$$c External functions
1537 c$$$ external ran_number
1538 c$$$ double precision ran_number
1540 c$$$c Input arguments
1542 c$$$ double precision angle ! In degrees
1544 c$$$c Local variables
1546 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1550 c$$$ rad_ang=angle*deg2rad
1553 c$$$ do while (length.lt.0.01)
1554 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1555 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1556 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1557 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1558 c$$$ + rand_v(3)*rand_v(3)
1559 c$$$ length=sqrt(length)
1560 c$$$ rand_v(1)=rand_v(1)/length
1561 c$$$ rand_v(2)=rand_v(2)/length
1562 c$$$ rand_v(3)=rand_v(3)/length
1563 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1564 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1565 c$$$ length=1.0D0-cost*cost
1566 c$$$ if (length.lt.0.0D0) length=0.0D0
1567 c$$$ length=sqrt(length)
1568 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1569 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1570 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1572 c$$$ rand_v(1)=rand_v(1)/length
1573 c$$$ rand_v(2)=rand_v(2)/length
1574 c$$$ rand_v(3)=rand_v(3)/length
1576 c$$$ cost=dcos(rad_ang)
1577 c$$$ sint=dsin(rad_ang)
1578 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1579 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1580 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1581 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1582 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1583 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1584 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1585 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1586 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1588 c$$$ call chainbuild_cart
1593 c$$$c----------------------------------------------------------------------------
1595 c$$$ subroutine ss_relax3(i_in,j_in)
1599 c$$$ include 'DIMENSIONS'
1600 c$$$ include 'COMMON.VAR'
1601 c$$$ include 'COMMON.CHAIN'
1602 c$$$ include 'COMMON.IOUNITS'
1603 c$$$ include 'COMMON.INTERACT'
1605 c$$$c External functions
1606 c$$$ external ran_number
1607 c$$$ double precision ran_number
1609 c$$$c Input arguments
1610 c$$$ integer i_in,j_in
1612 c$$$c Local variables
1613 c$$$ double precision energy_sc(0:n_ene),etot
1614 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1615 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1616 c$$$ integer n,i_pert,i
1617 c$$$ logical notdone
1626 c$$$ mask_side(i_in)=1
1627 c$$$ mask_side(j_in)=1
1629 c$$$ call etotal_sc(energy_sc)
1630 c$$$ etot=energy_sc(0)
1631 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1632 c$$$c + energy_sc(1),energy_sc(12)
1636 c$$$ do while (notdone)
1637 c$$$ if (mod(n,2).eq.0) then
1645 c$$$ org_dc(i)=dc(i,i_pert+nres)
1646 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1647 c$$$ org_c(i)=c(i,i_pert+nres)
1649 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1650 c$$$ call perturb_side_chain(i_pert,ang_pert)
1651 c$$$ call etotal_sc(energy_sc)
1652 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1653 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1654 c$$$ if (rand_fact.lt.exp_fact) then
1655 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1656 c$$$c + energy_sc(1),energy_sc(12)
1657 c$$$ etot=energy_sc(0)
1659 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1660 c$$$c + energy_sc(1),energy_sc(12)
1662 c$$$ dc(i,i_pert+nres)=org_dc(i)
1663 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1664 c$$$ c(i,i_pert+nres)=org_c(i)
1668 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1676 c$$$c----------------------------------------------------------------------------
1678 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1680 c$$$ include 'DIMENSIONS'
1682 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1683 c$$$*********************************************************************
1684 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1685 c$$$* the calling subprogram. *
1686 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1687 c$$$* calculated in the usual pythagorean way. *
1688 c$$$* absolute convergence occurs when the function is within v(31) of *
1689 c$$$* zero. unless you know the minimum value in advance, abs convg *
1690 c$$$* is probably not useful. *
1691 c$$$* relative convergence is when the model predicts that the function *
1692 c$$$* will decrease by less than v(32)*abs(fun). *
1693 c$$$*********************************************************************
1694 c$$$ include 'COMMON.IOUNITS'
1695 c$$$ include 'COMMON.VAR'
1696 c$$$ include 'COMMON.GEO'
1697 c$$$ include 'COMMON.MINIM'
1698 c$$$ include 'COMMON.CHAIN'
1700 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1701 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1702 c$$$ + orig_ss_dist(maxres2,maxres2)
1704 c$$$ double precision etot
1705 c$$$ integer iretcode,nfun,i_in,j_in
1708 c$$$ double precision dist
1709 c$$$ external ss_func,fdum
1710 c$$$ double precision ss_func,fdum
1712 c$$$ integer iv(liv),uiparm(2)
1713 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1717 c$$$ call deflt(2,iv,liv,lv,v)
1718 c$$$* 12 means fresh start, dont call deflt
1720 c$$$* max num of fun calls
1721 c$$$ if (maxfun.eq.0) maxfun=500
1723 c$$$* max num of iterations
1724 c$$$ if (maxmin.eq.0) maxmin=1000
1726 c$$$* controls output
1728 c$$$* selects output unit
1731 c$$$* 1 means to print out result
1733 c$$$* 1 means to print out summary stats
1735 c$$$* 1 means to print initial x and d
1737 c$$$* min val for v(radfac) default is 0.1
1739 c$$$* max val for v(radfac) default is 4.0
1742 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1743 c$$$* the sumsl default is 0.1
1745 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1746 c$$$* the sumsl default is 100*machep
1747 c$$$ v(34)=v(34)/100.0D0
1748 c$$$* absolute convergence
1749 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1752 c$$$* relative convergence
1753 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1756 c$$$* controls initial step size
1758 c$$$* large vals of d correspond to small components of step
1765 c$$$ orig_ss_dc(j,i)=dc(j,i)
1768 c$$$ call geom_to_var(nvar,orig_ss_var)
1772 c$$$ orig_ss_dist(j,i)=dist(j,i)
1773 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1774 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1775 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1787 c$$$ if (ialph(i,1).gt.0) then
1790 c$$$ x(k)=dc(j,i+nres)
1797 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1800 c$$$ nfun=iv(6)+iv(30)
1810 c$$$ if (ialph(i,1).gt.0) then
1813 c$$$ dc(j,i+nres)=x(k)
1817 c$$$ call chainbuild_cart
1822 c$$$C-----------------------------------------------------------------------------
1824 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1826 c$$$ include 'DIMENSIONS'
1827 c$$$ include 'COMMON.DERIV'
1828 c$$$ include 'COMMON.IOUNITS'
1829 c$$$ include 'COMMON.VAR'
1830 c$$$ include 'COMMON.CHAIN'
1831 c$$$ include 'COMMON.INTERACT'
1832 c$$$ include 'COMMON.SBRIDGE'
1834 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1835 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1836 c$$$ + orig_ss_dist(maxres2,maxres2)
1839 c$$$ double precision x(maxres6)
1841 c$$$ double precision f
1842 c$$$ integer uiparm(2)
1843 c$$$ real*8 urparm(1)
1844 c$$$ external ufparm
1845 c$$$ double precision ufparm
1848 c$$$ double precision dist
1850 c$$$ integer i,j,k,ss_i,ss_j
1851 c$$$ double precision tempf,var(maxvar)
1866 c$$$ if (ialph(i,1).gt.0) then
1869 c$$$ dc(j,i+nres)=x(k)
1873 c$$$ call chainbuild_cart
1875 c$$$ call geom_to_var(nvar,var)
1877 c$$$c Constraints on all angles
1879 c$$$ tempf=var(i)-orig_ss_var(i)
1880 c$$$ f=f+tempf*tempf
1883 c$$$c Constraints on all distances
1885 c$$$ if (i.gt.1) then
1886 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1887 c$$$ f=f+tempf*tempf
1890 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1891 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1892 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1893 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1894 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1895 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1896 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1897 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1901 c$$$c Constraints for the relevant CYS-CYS
1902 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1903 c$$$ f=f+tempf*tempf
1904 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
1906 c$$$c$$$ if (nf.ne.nfl) then
1907 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1908 c$$$c$$$ + f,dist(5+nres,14+nres)
1916 c$$$C-----------------------------------------------------------------------------