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(0: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 c if (.not.found.and.fg_rank.eq.0)
628 c & write(iout,'(a15,f12.2,f8.1,2i5)')
629 c & "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 c if (.not.found.and.fg_rank.eq.0)
643 c & write(iout,'(a15,f12.2,f8.1,2i5)')
644 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
659 C-----------------------------------------------------------------------------
660 C-----------------------------------------------------------------------------
661 C-----------------------------------------------------------------------------
663 c$$$c-----------------------------------------------------------------------------
665 c$$$ subroutine ss_relax(i_in,j_in)
669 c$$$ include 'DIMENSIONS'
670 c$$$ include 'COMMON.VAR'
671 c$$$ include 'COMMON.CHAIN'
672 c$$$ include 'COMMON.IOUNITS'
673 c$$$ include 'COMMON.INTERACT'
675 c$$$c Input arguments
676 c$$$ integer i_in,j_in
678 c$$$c Local variables
679 c$$$ integer i,iretcode,nfun_sc
681 c$$$ double precision var(maxvar),e_sc,etot
688 c$$$ mask_side(i_in)=1
689 c$$$ mask_side(j_in)=1
691 c$$$c Minimize the two selected side-chains
692 c$$$ call overlap_sc(scfail) ! Better not fail!
693 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
700 c$$$c-------------------------------------------------------------
702 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
703 c$$$c Minimize side-chains only, starting from geom but without modifying
705 c$$$c If mask_r is already set, only the selected side-chains are minimized,
706 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
710 c$$$ include 'DIMENSIONS'
711 c$$$ include 'COMMON.IOUNITS'
712 c$$$ include 'COMMON.VAR'
713 c$$$ include 'COMMON.CHAIN'
714 c$$$ include 'COMMON.GEO'
715 c$$$ include 'COMMON.MINIM'
717 c$$$ common /srutu/ icall
719 c$$$c Output arguments
720 c$$$ double precision etot_sc
721 c$$$ integer iretcode,nfun
723 c$$$c External functions/subroutines
724 c$$$ external func_sc,grad_sc,fdum
726 c$$$c Local variables
728 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
730 c$$$ double precision rdum(1)
731 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
733 c$$$ integer i,nvar_restr
736 c$$$cmc start_minim=.true.
737 c$$$ call deflt(2,iv,liv,lv,v)
738 c$$$* 12 means fresh start, dont call deflt
740 c$$$* max num of fun calls
741 c$$$ if (maxfun.eq.0) maxfun=500
743 c$$$* max num of iterations
744 c$$$ if (maxmin.eq.0) maxmin=1000
746 c$$$* controls output
748 c$$$* selects output unit
750 c$$$c iv(21)=iout ! DEBUG
751 c$$$c iv(21)=8 ! DEBUG
752 c$$$* 1 means to print out result
754 c$$$c iv(22)=1 ! DEBUG
755 c$$$* 1 means to print out summary stats
757 c$$$c iv(23)=1 ! DEBUG
758 c$$$* 1 means to print initial x and d
760 c$$$c iv(24)=1 ! DEBUG
761 c$$$* min val for v(radfac) default is 0.1
763 c$$$* max val for v(radfac) default is 4.0
766 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
767 c$$$* the sumsl default is 0.1
769 c$$$* false conv if (act fnctn decrease) .lt. v(34)
770 c$$$* the sumsl default is 100*machep
771 c$$$ v(34)=v(34)/100.0D0
772 c$$$* absolute convergence
773 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
775 c$$$* relative convergence
776 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
778 c$$$* controls initial step size
780 c$$$* large vals of d correspond to small components of step
784 c$$$ do i=nphi+1,nvar
788 c$$$ call geom_to_var(nvar,x)
789 c$$$ IF (mask_r) THEN
790 c$$$ do i=1,nres ! Just in case...
794 c$$$ call x2xx(x,xx,nvar_restr)
795 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
796 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
799 c$$$c When minimizing ALL side-chains, etotal_sc is a little
800 c$$$c faster if we don't set mask_r
806 c$$$ call x2xx(x,xx,nvar_restr)
807 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
808 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
811 c$$$ call var_to_geom(nvar,x)
812 c$$$ call chainbuild_sc
819 c$$$C--------------------------------------------------------------------------
821 c$$$ subroutine chainbuild_sc
823 c$$$ include 'DIMENSIONS'
824 c$$$ include 'COMMON.VAR'
825 c$$$ include 'COMMON.INTERACT'
827 c$$$c Local variables
832 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
833 c$$$ call locate_side_chain(i)
840 c$$$C--------------------------------------------------------------------------
842 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
846 c$$$ include 'DIMENSIONS'
847 c$$$ include 'COMMON.DERIV'
848 c$$$ include 'COMMON.VAR'
849 c$$$ include 'COMMON.MINIM'
850 c$$$ include 'COMMON.IOUNITS'
852 c$$$c Input arguments
854 c$$$ double precision x(maxvar)
855 c$$$ double precision ufparm
858 c$$$c Input/Output arguments
860 c$$$ integer uiparm(1)
861 c$$$ double precision urparm(1)
863 c$$$c Output arguments
864 c$$$ double precision f
866 c$$$c Local variables
867 c$$$ double precision energia(0:n_ene)
869 c$$$c Variables used to intercept NaNs
870 c$$$ double precision x_sum
879 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
882 c$$$ x_sum=x_sum+x(i_NAN)
884 c$$$c Calculate the energy only if the coordinates are ok
885 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
886 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
892 c$$$ call var_to_geom_restr(n,x)
894 c$$$ call chainbuild_sc
895 c$$$ call etotal_sc(energia(0))
897 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
906 c$$$c-------------------------------------------------------
908 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
912 c$$$ include 'DIMENSIONS'
913 c$$$ include 'COMMON.CHAIN'
914 c$$$ include 'COMMON.DERIV'
915 c$$$ include 'COMMON.VAR'
916 c$$$ include 'COMMON.INTERACT'
917 c$$$ include 'COMMON.MINIM'
919 c$$$c Input arguments
921 c$$$ double precision x(maxvar)
922 c$$$ double precision ufparm
925 c$$$c Input/Output arguments
927 c$$$ integer uiparm(1)
928 c$$$ double precision urparm(1)
930 c$$$c Output arguments
931 c$$$ double precision g(maxvar)
933 c$$$c Local variables
934 c$$$ double precision f,gphii,gthetai,galphai,gomegai
935 c$$$ integer ig,ind,i,j,k,igall,ij
939 c$$$ if (nf-nfl+1) 20,30,40
940 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
941 c$$$c write (iout,*) 'grad 20'
942 c$$$ if (nf.eq.0) return
944 c$$$ 30 call var_to_geom_restr(n,x)
945 c$$$ call chainbuild_sc
947 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
951 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
957 c$$$ IF (mask_phi(i+2).eq.1) THEN
962 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
963 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
969 c$$$ ind=ind+nres-1-i
976 c$$$ IF (mask_theta(i+2).eq.1) THEN
982 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
983 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
988 c$$$ ind=ind+nres-1-i
993 c$$$ if (itype(i).ne.10) then
994 c$$$ IF (mask_side(i).eq.1) THEN
998 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1007 c$$$ if (itype(i).ne.10) then
1008 c$$$ IF (mask_side(i).eq.1) THEN
1012 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1020 c$$$C Add the components corresponding to local energy terms.
1027 c$$$ if (mask_phi(i).eq.1) then
1029 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1035 c$$$ if (mask_theta(i).eq.1) then
1037 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1043 c$$$ if (itype(i).ne.10) then
1045 c$$$ if (mask_side(i).eq.1) then
1047 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1054 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1060 c$$$C-----------------------------------------------------------------------------
1062 c$$$ subroutine etotal_sc(energy_sc)
1066 c$$$ include 'DIMENSIONS'
1067 c$$$ include 'COMMON.VAR'
1068 c$$$ include 'COMMON.INTERACT'
1069 c$$$ include 'COMMON.DERIV'
1070 c$$$ include 'COMMON.FFIELD'
1072 c$$$c Output arguments
1073 c$$$ double precision energy_sc(0:n_ene)
1075 c$$$c Local variables
1076 c$$$ double precision evdw,escloc
1081 c$$$ energy_sc(i)=0.0D0
1084 c$$$ if (mask_r) then
1085 c$$$ call egb_sc(evdw)
1086 c$$$ call esc_sc(escloc)
1089 c$$$ call esc(escloc)
1092 c$$$ if (evdw.eq.1.0D20) then
1093 c$$$ energy_sc(0)=evdw
1095 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1097 c$$$ energy_sc(1)=evdw
1098 c$$$ energy_sc(12)=escloc
1101 c$$$C Sum up the components of the Cartesian gradient.
1105 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1112 c$$$C-----------------------------------------------------------------------------
1114 c$$$ subroutine egb_sc(evdw)
1116 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1117 c$$$C assuming the Gay-Berne potential of interaction.
1119 c$$$ implicit real*8 (a-h,o-z)
1120 c$$$ include 'DIMENSIONS'
1121 c$$$ include 'COMMON.GEO'
1122 c$$$ include 'COMMON.VAR'
1123 c$$$ include 'COMMON.LOCAL'
1124 c$$$ include 'COMMON.CHAIN'
1125 c$$$ include 'COMMON.DERIV'
1126 c$$$ include 'COMMON.NAMES'
1127 c$$$ include 'COMMON.INTERACT'
1128 c$$$ include 'COMMON.IOUNITS'
1129 c$$$ include 'COMMON.CALC'
1130 c$$$ include 'COMMON.CONTROL'
1133 c$$$ energy_dec=.false.
1134 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1137 c$$$c if (icall.eq.0) lprn=.false.
1139 c$$$ do i=iatsc_s,iatsc_e
1141 c$$$ itypi1=itype(i+1)
1145 c$$$ dxi=dc_norm(1,nres+i)
1146 c$$$ dyi=dc_norm(2,nres+i)
1147 c$$$ dzi=dc_norm(3,nres+i)
1148 c$$$c dsci_inv=dsc_inv(itypi)
1149 c$$$ dsci_inv=vbld_inv(i+nres)
1150 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1151 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1153 c$$$C Calculate SC interaction energy.
1155 c$$$ do iint=1,nint_gr(i)
1156 c$$$ do j=istart(i,iint),iend(i,iint)
1157 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1160 c$$$c dscj_inv=dsc_inv(itypj)
1161 c$$$ dscj_inv=vbld_inv(j+nres)
1162 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1163 c$$$c & 1.0d0/vbld(j+nres)
1164 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1165 c$$$ sig0ij=sigma(itypi,itypj)
1166 c$$$ chi1=chi(itypi,itypj)
1167 c$$$ chi2=chi(itypj,itypi)
1168 c$$$ chi12=chi1*chi2
1169 c$$$ chip1=chip(itypi)
1170 c$$$ chip2=chip(itypj)
1171 c$$$ chip12=chip1*chip2
1172 c$$$ alf1=alp(itypi)
1173 c$$$ alf2=alp(itypj)
1174 c$$$ alf12=0.5D0*(alf1+alf2)
1175 c$$$C For diagnostics only!!!
1185 c$$$ xj=c(1,nres+j)-xi
1186 c$$$ yj=c(2,nres+j)-yi
1187 c$$$ zj=c(3,nres+j)-zi
1188 c$$$ dxj=dc_norm(1,nres+j)
1189 c$$$ dyj=dc_norm(2,nres+j)
1190 c$$$ dzj=dc_norm(3,nres+j)
1191 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1192 c$$$c write (iout,*) "j",j," dc_norm",
1193 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1194 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1195 c$$$ rij=dsqrt(rrij)
1196 c$$$C Calculate angle-dependent terms of energy and contributions to their
1198 c$$$ call sc_angular
1199 c$$$ sigsq=1.0D0/sigsq
1200 c$$$ sig=sig0ij*dsqrt(sigsq)
1201 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1202 c$$$c for diagnostics; uncomment
1203 c$$$c rij_shift=1.2*sig0ij
1204 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1205 c$$$ if (rij_shift.le.0.0D0) then
1207 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1208 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1209 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1212 c$$$ sigder=-sig*sigsq
1213 c$$$c---------------------------------------------------------------
1214 c$$$ rij_shift=1.0D0/rij_shift
1215 c$$$ fac=rij_shift**expon
1216 c$$$ e1=fac*fac*aa(itypi,itypj)
1217 c$$$ e2=fac*bb(itypi,itypj)
1218 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1219 c$$$ eps2der=evdwij*eps3rt
1220 c$$$ eps3der=evdwij*eps2rt
1221 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1222 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1223 c$$$ evdwij=evdwij*eps2rt*eps3rt
1224 c$$$ evdw=evdw+evdwij
1226 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1227 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1228 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1229 c$$$ & restyp(itypi),i,restyp(itypj),j,
1230 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1231 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1232 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1236 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1237 c$$$ & 'evdw',i,j,evdwij
1239 c$$$C Calculate gradient components.
1240 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1241 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1242 c$$$ sigder=fac*sigder
1245 c$$$C Calculate the radial part of the gradient
1249 c$$$C Calculate angular part of the gradient.
1255 c$$$ energy_dec=.false.
1259 c$$$c-----------------------------------------------------------------------------
1261 c$$$ subroutine esc_sc(escloc)
1262 c$$$C Calculate the local energy of a side chain and its derivatives in the
1263 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1264 c$$$C ALPHA and OMEGA.
1265 c$$$ implicit real*8 (a-h,o-z)
1266 c$$$ include 'DIMENSIONS'
1267 c$$$ include 'COMMON.GEO'
1268 c$$$ include 'COMMON.LOCAL'
1269 c$$$ include 'COMMON.VAR'
1270 c$$$ include 'COMMON.INTERACT'
1271 c$$$ include 'COMMON.DERIV'
1272 c$$$ include 'COMMON.CHAIN'
1273 c$$$ include 'COMMON.IOUNITS'
1274 c$$$ include 'COMMON.NAMES'
1275 c$$$ include 'COMMON.FFIELD'
1276 c$$$ include 'COMMON.CONTROL'
1277 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1278 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1279 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1280 c$$$ delta=0.02d0*pi
1282 c$$$c write (iout,'(a)') 'ESC'
1283 c$$$ do i=loc_start,loc_end
1284 c$$$ IF (mask_side(i).eq.1) THEN
1286 c$$$ if (it.eq.10) goto 1
1287 c$$$ nlobit=nlob(it)
1288 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1289 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1290 c$$$ theti=theta(i+1)-pipol
1291 c$$$ x(1)=dtan(theti)
1295 c$$$ if (x(2).gt.pi-delta) then
1297 c$$$ xtemp(2)=pi-delta
1299 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1301 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1302 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1303 c$$$ & escloci,dersc(2))
1304 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1305 c$$$ & ddersc0(1),dersc(1))
1306 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1307 c$$$ & ddersc0(3),dersc(3))
1308 c$$$ xtemp(2)=pi-delta
1309 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1311 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1312 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1313 c$$$ & dersc0(2),esclocbi,dersc02)
1314 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1315 c$$$ & dersc12,dersc01)
1316 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1317 c$$$ dersc0(1)=dersc01
1318 c$$$ dersc0(2)=dersc02
1319 c$$$ dersc0(3)=0.0d0
1321 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1323 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1324 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1325 c$$$c & esclocbi,ss,ssd
1326 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1327 c$$$c escloci=esclocbi
1328 c$$$c write (iout,*) escloci
1329 c$$$ else if (x(2).lt.delta) then
1333 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1335 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1336 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1337 c$$$ & escloci,dersc(2))
1338 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1339 c$$$ & ddersc0(1),dersc(1))
1340 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1341 c$$$ & ddersc0(3),dersc(3))
1343 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1345 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1346 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1347 c$$$ & dersc0(2),esclocbi,dersc02)
1348 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1349 c$$$ & dersc12,dersc01)
1350 c$$$ dersc0(1)=dersc01
1351 c$$$ dersc0(2)=dersc02
1352 c$$$ dersc0(3)=0.0d0
1353 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1355 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1357 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1358 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1359 c$$$c & esclocbi,ss,ssd
1360 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1361 c$$$c write (iout,*) escloci
1363 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1366 c$$$ escloc=escloc+escloci
1367 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1368 c$$$ & 'escloc',i,escloci
1369 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1371 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1372 c$$$ & wscloc*dersc(1)
1373 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1374 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1381 c$$$C-----------------------------------------------------------------------------
1383 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1385 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1386 c$$$C assuming the Gay-Berne potential of interaction.
1388 c$$$ implicit real*8 (a-h,o-z)
1389 c$$$ include 'DIMENSIONS'
1390 c$$$ include 'COMMON.GEO'
1391 c$$$ include 'COMMON.VAR'
1392 c$$$ include 'COMMON.LOCAL'
1393 c$$$ include 'COMMON.CHAIN'
1394 c$$$ include 'COMMON.DERIV'
1395 c$$$ include 'COMMON.NAMES'
1396 c$$$ include 'COMMON.INTERACT'
1397 c$$$ include 'COMMON.IOUNITS'
1398 c$$$ include 'COMMON.CALC'
1399 c$$$ include 'COMMON.CONTROL'
1402 c$$$ energy_dec=.false.
1403 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1407 c$$$c$$$ do i=iatsc_s,iatsc_e
1410 c$$$ itypi1=itype(i+1)
1414 c$$$ dxi=dc_norm(1,nres+i)
1415 c$$$ dyi=dc_norm(2,nres+i)
1416 c$$$ dzi=dc_norm(3,nres+i)
1417 c$$$c dsci_inv=dsc_inv(itypi)
1418 c$$$ dsci_inv=vbld_inv(i+nres)
1419 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1420 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1422 c$$$C Calculate SC interaction energy.
1424 c$$$c$$$ do iint=1,nint_gr(i)
1425 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1429 c$$$c dscj_inv=dsc_inv(itypj)
1430 c$$$ dscj_inv=vbld_inv(j+nres)
1431 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1432 c$$$c & 1.0d0/vbld(j+nres)
1433 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1434 c$$$ sig0ij=sigma(itypi,itypj)
1435 c$$$ chi1=chi(itypi,itypj)
1436 c$$$ chi2=chi(itypj,itypi)
1437 c$$$ chi12=chi1*chi2
1438 c$$$ chip1=chip(itypi)
1439 c$$$ chip2=chip(itypj)
1440 c$$$ chip12=chip1*chip2
1441 c$$$ alf1=alp(itypi)
1442 c$$$ alf2=alp(itypj)
1443 c$$$ alf12=0.5D0*(alf1+alf2)
1444 c$$$C For diagnostics only!!!
1454 c$$$ xj=c(1,nres+j)-xi
1455 c$$$ yj=c(2,nres+j)-yi
1456 c$$$ zj=c(3,nres+j)-zi
1457 c$$$ dxj=dc_norm(1,nres+j)
1458 c$$$ dyj=dc_norm(2,nres+j)
1459 c$$$ dzj=dc_norm(3,nres+j)
1460 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1461 c$$$c write (iout,*) "j",j," dc_norm",
1462 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1463 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1464 c$$$ rij=dsqrt(rrij)
1465 c$$$C Calculate angle-dependent terms of energy and contributions to their
1467 c$$$ call sc_angular
1468 c$$$ sigsq=1.0D0/sigsq
1469 c$$$ sig=sig0ij*dsqrt(sigsq)
1470 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1471 c$$$c for diagnostics; uncomment
1472 c$$$c rij_shift=1.2*sig0ij
1473 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1474 c$$$ if (rij_shift.le.0.0D0) then
1476 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1477 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1478 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1481 c$$$ sigder=-sig*sigsq
1482 c$$$c---------------------------------------------------------------
1483 c$$$ rij_shift=1.0D0/rij_shift
1484 c$$$ fac=rij_shift**expon
1485 c$$$ e1=fac*fac*aa(itypi,itypj)
1486 c$$$ e2=fac*bb(itypi,itypj)
1487 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1488 c$$$ eps2der=evdwij*eps3rt
1489 c$$$ eps3der=evdwij*eps2rt
1490 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1491 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1492 c$$$ evdwij=evdwij*eps2rt*eps3rt
1493 c$$$ evdw=evdw+evdwij
1495 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1496 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1497 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1498 c$$$ & restyp(itypi),i,restyp(itypj),j,
1499 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1500 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1501 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1505 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1506 c$$$ & 'evdw',i,j,evdwij
1508 c$$$C Calculate gradient components.
1509 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1510 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1511 c$$$ sigder=fac*sigder
1514 c$$$C Calculate the radial part of the gradient
1518 c$$$C Calculate angular part of the gradient.
1521 c$$$c$$$ enddo ! iint
1523 c$$$ energy_dec=.false.
1527 c$$$C-----------------------------------------------------------------------------
1529 c$$$ subroutine perturb_side_chain(i,angle)
1533 c$$$ include 'DIMENSIONS'
1534 c$$$ include 'COMMON.CHAIN'
1535 c$$$ include 'COMMON.GEO'
1536 c$$$ include 'COMMON.VAR'
1537 c$$$ include 'COMMON.LOCAL'
1538 c$$$ include 'COMMON.IOUNITS'
1540 c$$$c External functions
1541 c$$$ external ran_number
1542 c$$$ double precision ran_number
1544 c$$$c Input arguments
1546 c$$$ double precision angle ! In degrees
1548 c$$$c Local variables
1550 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1554 c$$$ rad_ang=angle*deg2rad
1557 c$$$ do while (length.lt.0.01)
1558 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1559 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1560 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1561 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1562 c$$$ + rand_v(3)*rand_v(3)
1563 c$$$ length=sqrt(length)
1564 c$$$ rand_v(1)=rand_v(1)/length
1565 c$$$ rand_v(2)=rand_v(2)/length
1566 c$$$ rand_v(3)=rand_v(3)/length
1567 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1568 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1569 c$$$ length=1.0D0-cost*cost
1570 c$$$ if (length.lt.0.0D0) length=0.0D0
1571 c$$$ length=sqrt(length)
1572 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1573 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1574 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1576 c$$$ rand_v(1)=rand_v(1)/length
1577 c$$$ rand_v(2)=rand_v(2)/length
1578 c$$$ rand_v(3)=rand_v(3)/length
1580 c$$$ cost=dcos(rad_ang)
1581 c$$$ sint=dsin(rad_ang)
1582 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1583 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1584 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1585 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1586 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1587 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1588 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1589 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1590 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1592 c$$$ call chainbuild_cart
1597 c$$$c----------------------------------------------------------------------------
1599 c$$$ subroutine ss_relax3(i_in,j_in)
1603 c$$$ include 'DIMENSIONS'
1604 c$$$ include 'COMMON.VAR'
1605 c$$$ include 'COMMON.CHAIN'
1606 c$$$ include 'COMMON.IOUNITS'
1607 c$$$ include 'COMMON.INTERACT'
1609 c$$$c External functions
1610 c$$$ external ran_number
1611 c$$$ double precision ran_number
1613 c$$$c Input arguments
1614 c$$$ integer i_in,j_in
1616 c$$$c Local variables
1617 c$$$ double precision energy_sc(0:n_ene),etot
1618 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1619 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1620 c$$$ integer n,i_pert,i
1621 c$$$ logical notdone
1630 c$$$ mask_side(i_in)=1
1631 c$$$ mask_side(j_in)=1
1633 c$$$ call etotal_sc(energy_sc)
1634 c$$$ etot=energy_sc(0)
1635 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1636 c$$$c + energy_sc(1),energy_sc(12)
1640 c$$$ do while (notdone)
1641 c$$$ if (mod(n,2).eq.0) then
1649 c$$$ org_dc(i)=dc(i,i_pert+nres)
1650 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1651 c$$$ org_c(i)=c(i,i_pert+nres)
1653 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1654 c$$$ call perturb_side_chain(i_pert,ang_pert)
1655 c$$$ call etotal_sc(energy_sc)
1656 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1657 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1658 c$$$ if (rand_fact.lt.exp_fact) then
1659 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1660 c$$$c + energy_sc(1),energy_sc(12)
1661 c$$$ etot=energy_sc(0)
1663 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1664 c$$$c + energy_sc(1),energy_sc(12)
1666 c$$$ dc(i,i_pert+nres)=org_dc(i)
1667 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1668 c$$$ c(i,i_pert+nres)=org_c(i)
1672 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1680 c$$$c----------------------------------------------------------------------------
1682 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1684 c$$$ include 'DIMENSIONS'
1686 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1687 c$$$*********************************************************************
1688 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1689 c$$$* the calling subprogram. *
1690 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1691 c$$$* calculated in the usual pythagorean way. *
1692 c$$$* absolute convergence occurs when the function is within v(31) of *
1693 c$$$* zero. unless you know the minimum value in advance, abs convg *
1694 c$$$* is probably not useful. *
1695 c$$$* relative convergence is when the model predicts that the function *
1696 c$$$* will decrease by less than v(32)*abs(fun). *
1697 c$$$*********************************************************************
1698 c$$$ include 'COMMON.IOUNITS'
1699 c$$$ include 'COMMON.VAR'
1700 c$$$ include 'COMMON.GEO'
1701 c$$$ include 'COMMON.MINIM'
1702 c$$$ include 'COMMON.CHAIN'
1704 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1705 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1706 c$$$ + orig_ss_dist(maxres2,maxres2)
1708 c$$$ double precision etot
1709 c$$$ integer iretcode,nfun,i_in,j_in
1712 c$$$ double precision dist
1713 c$$$ external ss_func,fdum
1714 c$$$ double precision ss_func,fdum
1716 c$$$ integer iv(liv),uiparm(2)
1717 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1721 c$$$ call deflt(2,iv,liv,lv,v)
1722 c$$$* 12 means fresh start, dont call deflt
1724 c$$$* max num of fun calls
1725 c$$$ if (maxfun.eq.0) maxfun=500
1727 c$$$* max num of iterations
1728 c$$$ if (maxmin.eq.0) maxmin=1000
1730 c$$$* controls output
1732 c$$$* selects output unit
1735 c$$$* 1 means to print out result
1737 c$$$* 1 means to print out summary stats
1739 c$$$* 1 means to print initial x and d
1741 c$$$* min val for v(radfac) default is 0.1
1743 c$$$* max val for v(radfac) default is 4.0
1746 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1747 c$$$* the sumsl default is 0.1
1749 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1750 c$$$* the sumsl default is 100*machep
1751 c$$$ v(34)=v(34)/100.0D0
1752 c$$$* absolute convergence
1753 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1756 c$$$* relative convergence
1757 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1760 c$$$* controls initial step size
1762 c$$$* large vals of d correspond to small components of step
1769 c$$$ orig_ss_dc(j,i)=dc(j,i)
1772 c$$$ call geom_to_var(nvar,orig_ss_var)
1776 c$$$ orig_ss_dist(j,i)=dist(j,i)
1777 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1778 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1779 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1791 c$$$ if (ialph(i,1).gt.0) then
1794 c$$$ x(k)=dc(j,i+nres)
1801 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1804 c$$$ nfun=iv(6)+iv(30)
1814 c$$$ if (ialph(i,1).gt.0) then
1817 c$$$ dc(j,i+nres)=x(k)
1821 c$$$ call chainbuild_cart
1826 c$$$C-----------------------------------------------------------------------------
1828 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1830 c$$$ include 'DIMENSIONS'
1831 c$$$ include 'COMMON.DERIV'
1832 c$$$ include 'COMMON.IOUNITS'
1833 c$$$ include 'COMMON.VAR'
1834 c$$$ include 'COMMON.CHAIN'
1835 c$$$ include 'COMMON.INTERACT'
1836 c$$$ include 'COMMON.SBRIDGE'
1838 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1839 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1840 c$$$ + orig_ss_dist(maxres2,maxres2)
1843 c$$$ double precision x(maxres6)
1845 c$$$ double precision f
1846 c$$$ integer uiparm(2)
1847 c$$$ real*8 urparm(1)
1848 c$$$ external ufparm
1849 c$$$ double precision ufparm
1852 c$$$ double precision dist
1854 c$$$ integer i,j,k,ss_i,ss_j
1855 c$$$ double precision tempf,var(maxvar)
1870 c$$$ if (ialph(i,1).gt.0) then
1873 c$$$ dc(j,i+nres)=x(k)
1877 c$$$ call chainbuild_cart
1879 c$$$ call geom_to_var(nvar,var)
1881 c$$$c Constraints on all angles
1883 c$$$ tempf=var(i)-orig_ss_var(i)
1884 c$$$ f=f+tempf*tempf
1887 c$$$c Constraints on all distances
1889 c$$$ if (i.gt.1) then
1890 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1891 c$$$ f=f+tempf*tempf
1894 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1895 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1896 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1897 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1898 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1899 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1900 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1901 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1905 c$$$c Constraints for the relevant CYS-CYS
1906 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1907 c$$$ f=f+tempf*tempf
1908 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
1910 c$$$c$$$ if (nf.ne.nfl) then
1911 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1912 c$$$c$$$ + f,dist(5+nres,14+nres)
1920 c$$$C-----------------------------------------------------------------------------
1921 c$$$C-----------------------------------------------------------------------------
1922 subroutine triple_ssbond_ene(resi,resj,resk,eij)
1923 include 'DIMENSIONS'
1924 include 'COMMON.SBRIDGE'
1925 include 'COMMON.CHAIN'
1926 include 'COMMON.DERIV'
1927 include 'COMMON.LOCAL'
1928 include 'COMMON.INTERACT'
1929 include 'COMMON.VAR'
1930 include 'COMMON.IOUNITS'
1931 include 'COMMON.CALC'
1938 c External functions
1939 double precision h_base
1943 integer resi,resj,resk
1946 double precision eij,eij1,eij2,eij3
1950 c integer itypi,itypj,k,l
1951 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
1952 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
1953 double precision xik,yik,zik,xjk,yjk,zjk
1954 double precision sig0ij,ljd,sig,fac,e1,e2
1955 double precision dcosom1(3),dcosom2(3),ed
1956 double precision pom1,pom2
1957 double precision ljA,ljB,ljXs
1958 double precision d_ljB(1:3)
1959 double precision ssA,ssB,ssC,ssXs
1960 double precision ssxm,ljxm,ssm,ljm
1961 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
1963 if (dtriss.eq.0) return
1967 C write(iout,*) resi,resj,resk
1969 dxi=dc_norm(1,nres+i)
1970 dyi=dc_norm(2,nres+i)
1971 dzi=dc_norm(3,nres+i)
1972 dsci_inv=vbld_inv(i+nres)
1982 dxj=dc_norm(1,nres+j)
1983 dyj=dc_norm(2,nres+j)
1984 dzj=dc_norm(3,nres+j)
1985 dscj_inv=vbld_inv(j+nres)
1991 dxk=dc_norm(1,nres+k)
1992 dyk=dc_norm(2,nres+k)
1993 dzk=dc_norm(3,nres+k)
1994 dscj_inv=vbld_inv(k+nres)
2004 rrij=(xij*xij+yij*yij+zij*zij)
2005 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2006 rrik=(xik*xik+yik*yik+zik*zik)
2008 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2010 C there are three combination of distances for each trisulfide bonds
2011 C The first case the ith atom is the center
2012 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2013 C distance y is second distance the a,b,c,d are parameters derived for
2014 C this problem d parameter was set as a penalty currenlty set to 1.
2015 if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then
2018 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss)
2020 C second case jth atom is center
2021 if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then
2024 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss)
2026 C the third case kth atom is the center
2027 if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then
2030 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss)
2036 C write(iout,*)i,j,k,eij
2037 C The energy penalty calculated now time for the gradient part
2038 C derivative over rij
2039 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5)
2040 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)
2045 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2046 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2049 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2050 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2052 C now derivative over rik
2053 fac=-eij1**2/dtriss*
2054 &(-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5)
2055 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
2060 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2061 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2064 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2065 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2067 C now derivative over rjk
2068 fac=-eij2**2/dtriss*
2069 &(-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)-
2070 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
2075 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2076 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2079 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2080 gvdwc(l,k)=gvdwc(l,k)+gg(l)