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-----------------------------------------------------------------------------
1917 subroutine triple_ssbond_ene(resi,resj,resk,eij)
1918 include 'DIMENSIONS'
1919 include 'COMMON.SBRIDGE'
1920 include 'COMMON.CHAIN'
1921 include 'COMMON.DERIV'
1922 include 'COMMON.LOCAL'
1923 include 'COMMON.INTERACT'
1924 include 'COMMON.VAR'
1925 include 'COMMON.IOUNITS'
1926 include 'COMMON.CALC'
1933 c External functions
1934 double precision h_base
1938 integer resi,resj,resk
1941 double precision eij,eij1,eij2,eij3
1945 c integer itypi,itypj,k,l
1946 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
1947 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
1948 double precision xik,yik,zik,xjk,yjk,zjk
1949 double precision sig0ij,ljd,sig,fac,e1,e2
1950 double precision dcosom1(3),dcosom2(3),ed
1951 double precision pom1,pom2
1952 double precision ljA,ljB,ljXs
1953 double precision d_ljB(1:3)
1954 double precision ssA,ssB,ssC,ssXs
1955 double precision ssxm,ljxm,ssm,ljm
1956 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
1961 C write(iout,*) resi,resj,resk
1963 dxi=dc_norm(1,nres+i)
1964 dyi=dc_norm(2,nres+i)
1965 dzi=dc_norm(3,nres+i)
1966 dsci_inv=vbld_inv(i+nres)
1976 dxj=dc_norm(1,nres+j)
1977 dyj=dc_norm(2,nres+j)
1978 dzj=dc_norm(3,nres+j)
1979 dscj_inv=vbld_inv(j+nres)
1985 dxk=dc_norm(1,nres+k)
1986 dyk=dc_norm(2,nres+k)
1987 dzk=dc_norm(3,nres+k)
1988 dscj_inv=vbld_inv(k+nres)
1998 rrij=(xij*xij+yij*yij+zij*zij)
1999 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2000 rrik=(xik*xik+yik*yik+zik*zik)
2002 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2004 C there are three combination of distances for each trisulfide bonds
2005 C The first case the ith atom is the center
2006 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2007 C distance y is second distance the a,b,c,d are parameters derived for
2008 C this problem d parameter was set as a penalty currenlty set to 1.
2009 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2010 C second case jth atom is center
2011 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2012 C the third case kth atom is the center
2013 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2018 C write(iout,*)i,j,k,eij
2019 C The energy penalty calculated now time for the gradient part
2020 C derivative over rij
2021 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2022 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2027 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2028 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2031 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2032 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2034 C now derivative over rik
2035 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2036 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2041 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2042 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2045 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2046 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2048 C now derivative over rjk
2049 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2050 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2055 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2056 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2059 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2060 gvdwc(l,k)=gvdwc(l,k)+gg(l)