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'
101 C include 'COMMON.MD'
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)
156 C returning the ith atom to box
158 if (xi.lt.0) xi=xi+boxxsize
160 if (yi.lt.0) yi=yi+boxysize
162 if (zi.lt.0) zi=zi+boxzsize
163 if ((zi.gt.bordlipbot)
164 &.and.(zi.lt.bordliptop)) then
165 C the energy transfer exist
166 if (zi.lt.buflipbot) then
167 C what fraction I am in
169 & ((zi-bordlipbot)/lipbufthick)
170 C lipbufthick is thickenes of lipid buffore
171 sslipi=sscalelip(fracinbuf)
172 ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
173 elseif (zi.gt.bufliptop) then
174 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
175 sslipi=sscalelip(fracinbuf)
176 ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
188 C returning jth atom to box
190 if (xj.lt.0) xj=xj+boxxsize
192 if (yj.lt.0) yj=yj+boxysize
194 if (zj.lt.0) zj=zj+boxzsize
195 if ((zj.gt.bordlipbot)
196 &.and.(zj.lt.bordliptop)) then
197 C the energy transfer exist
198 if (zj.lt.buflipbot) then
199 C what fraction I am in
201 & ((zj-bordlipbot)/lipbufthick)
202 C lipbufthick is thickenes of lipid buffore
203 sslipj=sscalelip(fracinbuf)
204 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
205 elseif (zj.gt.bufliptop) then
206 fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
207 sslipj=sscalelip(fracinbuf)
208 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
217 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
218 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
219 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
220 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
223 C checking the distance
224 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
229 C finding the closest
233 xj=xj_safe+xshift*boxxsize
234 yj=yj_safe+yshift*boxysize
235 zj=zj_safe+zshift*boxzsize
236 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
237 if(dist_temp.lt.dist_init) then
247 if (subchap.eq.1) then
257 C xj=c(1,nres+j)-c(1,nres+i)
258 C yj=c(2,nres+j)-c(2,nres+i)
259 C zj=c(3,nres+j)-c(3,nres+i)
260 dxj=dc_norm(1,nres+j)
261 dyj=dc_norm(2,nres+j)
262 dzj=dc_norm(3,nres+j)
263 dscj_inv=vbld_inv(j+nres)
265 chi1=chi(itypi,itypj)
266 chi2=chi(itypj,itypi)
273 alf12=0.5D0*(alf1+alf2)
275 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
276 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
277 c The following are set in sc_angular
281 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
282 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
283 c om12=dxi*dxj+dyi*dyj+dzi*dzj
285 rij=1.0D0/rij ! Reset this so it makes sense
287 sig0ij=sigma(itypi,itypj)
288 sig=sig0ij*dsqrt(1.0D0/sigsq)
291 ljA=eps1*eps2rt**2*eps3rt**2
294 ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
299 deltat12=om2-om1+2.0d0
304 & +akth*(deltat1*deltat1+deltat2*deltat2)
305 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
306 ssxm=ssXs-0.5D0*ssB/ssA
309 c$$$c Some extra output
310 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
311 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
312 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
313 c$$$ if (ssx0.gt.0.0d0) then
314 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
318 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
319 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
320 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
322 c-------END TESTING CODE
325 c Stop and plot energy and derivative as a function of distance
327 ssm=ssC-0.25D0*ssB*ssB/ssA
328 ljm=-0.25D0*ljB*bb/aa
330 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
338 if (.not.checkstop) then
345 if (checkstop) rij=(ssxm-1.0d0)+
346 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
347 c-------END TESTING CODE
349 if (rij.gt.ljxm) then
352 fac=(1.0D0/ljd)**expon
355 eij=eps1*eps2rt*eps3rt*(e1+e2)
356 C write(iout,*) eij,'TU?1'
359 eij=eij*eps2rt*eps3rt
362 e1=e1*eps1*eps2rt**2*eps3rt**2
363 ed=-expon*(e1+eij)/ljd
365 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
366 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
367 eom12=eij*eps1_om12+eps2der*eps2rt_om12
368 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
369 else if (rij.lt.ssxm) then
372 eij=ssA*ssd*ssd+ssB*ssd+ssC
373 C write(iout,*) 'TU?2',ssc,ssd
374 ed=2*akcm*ssd+akct*deltat12
376 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
377 eom1=-2*akth*deltat1-pom1-om2*pom2
378 eom2= 2*akth*deltat2+pom1-om1*pom2
381 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
383 d_ssxm(1)=0.5D0*akct/ssA
387 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
388 d_ljxm(2)=d_ljxm(1)*sigsq_om2
389 d_ljxm(3)=d_ljxm(1)*sigsq_om12
390 d_ljxm(1)=d_ljxm(1)*sigsq_om1
392 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
395 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
399 ssm=ssC-0.25D0*ssB*ssB/ssA
400 d_ssm(1)=0.5D0*akct*ssB/ssA
401 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
402 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
404 f1=(rij-xm)/(ssxm-xm)
405 f2=(rij-ssxm)/(xm-ssxm)
409 C write(iout,*) eij,'TU?3'
410 delta_inv=1.0d0/(xm-ssxm)
411 deltasq_inv=delta_inv*delta_inv
413 fac1=deltasq_inv*fac*(xm-rij)
414 fac2=deltasq_inv*fac*(rij-ssxm)
415 ed=delta_inv*(Ht*hd2-ssm*hd1)
416 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
417 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
418 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
421 ljm=-0.25D0*ljB*bb/aa
422 d_ljm(1)=-0.5D0*bb/aa*ljB
423 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
424 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
426 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
427 f1=(rij-ljxm)/(xm-ljxm)
428 f2=(rij-xm)/(ljxm-xm)
432 C write(iout,*) 'TU?4',ssA
433 delta_inv=1.0d0/(ljxm-xm)
434 deltasq_inv=delta_inv*delta_inv
436 fac1=deltasq_inv*fac*(ljxm-rij)
437 fac2=deltasq_inv*fac*(rij-xm)
438 ed=delta_inv*(ljm*hd2-Ht*hd1)
439 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
440 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
441 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
443 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
445 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
451 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
452 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
453 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
455 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
456 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
457 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
458 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
461 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
463 c$$$ d_ljm(k)=ljm*d_ljB(k)
467 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
468 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
469 c$$$ d_ss(2)=akct*ssd
470 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
471 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
474 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
475 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
476 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
478 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
479 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
481 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
483 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
484 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
485 c$$$ h1=h_base(f1,hd1)
486 c$$$ h2=h_base(f2,hd2)
487 c$$$ eij=ss*h1+ljf*h2
488 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
489 c$$$ deltasq_inv=delta_inv*delta_inv
490 c$$$ fac=ljf*hd2-ss*hd1
491 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
492 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
493 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
494 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
495 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
496 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
497 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
499 c$$$ havebond=.false.
500 c$$$ if (ed.gt.0.0d0) havebond=.true.
501 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
504 write(iout,*) 'havebond',havebond
508 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
509 c write(iout,'(a15,f12.2,f8.1,2i5)')
510 c & "SSBOND_E_FORM",totT,t_bath,i,j
514 dyn_ssbond_ij(i,j)=eij
515 else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
516 dyn_ssbond_ij(i,j)=1.0d300
519 c write(iout,'(a15,f12.2,f8.1,2i5)')
520 c & "SSBOND_E_BREAK",totT,t_bath,i,j
527 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
528 & "CHECKSTOP",rij,eij,ed
533 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
540 c-------END TESTING CODE
543 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
544 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
547 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
550 gvdwx(k,i)=gvdwx(k,i)-gg(k)
551 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
552 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
553 gvdwx(k,j)=gvdwx(k,j)+gg(k)
554 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
555 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
559 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
564 gvdwc(l,i)=gvdwc(l,i)-gg(l)
565 gvdwc(l,j)=gvdwc(l,j)+gg(l)
571 C-----------------------------------------------------------------------------
573 double precision function h_base(x,deriv)
574 c A smooth function going 0->1 in range [0,1]
575 c It should NOT be called outside range [0,1], it will not work there.
582 double precision deriv
588 c Two parabolas put together. First derivative zero at extrema
589 c$$$ if (x.lt.0.5D0) then
590 c$$$ h_base=2.0D0*x*x
594 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
595 c$$$ deriv=4.0D0*deriv
598 c Third degree polynomial. First derivative zero at extrema
599 h_base=x*x*(3.0d0-2.0d0*x)
600 deriv=6.0d0*x*(1.0d0-x)
602 c Fifth degree polynomial. First and second derivatives zero at extrema
604 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
606 c$$$ deriv=deriv*deriv
607 c$$$ deriv=30.0d0*xsq*deriv
612 c----------------------------------------------------------------------------
614 subroutine dyn_set_nss
615 c Adjust nss and other relevant variables based on dyn_ssbond_ij
623 include 'COMMON.SBRIDGE'
624 include 'COMMON.CHAIN'
625 include 'COMMON.IOUNITS'
626 C include 'COMMON.SETUP'
629 C include 'COMMON.MD'
634 double precision emin
636 integer diff,allflag(maxdim),allnss,
637 & allihpb(maxdim),alljhpb(maxdim),
638 & newnss,newihpb(maxdim),newjhpb(maxdim)
640 integer i_newnss(1024),displ(0:1024)
641 integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
646 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
655 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
659 if (allflag(i).eq.0 .and.
660 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
661 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
665 if (emin.lt.1.0d300) then
668 if (allflag(i).eq.0 .and.
669 & (allihpb(i).eq.allihpb(imin) .or.
670 & alljhpb(i).eq.allihpb(imin) .or.
671 & allihpb(i).eq.alljhpb(imin) .or.
672 & alljhpb(i).eq.alljhpb(imin))) then
679 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
683 if (allflag(i).eq.1) then
685 newihpb(newnss)=allihpb(i)
686 newjhpb(newnss)=alljhpb(i)
691 if (nfgtasks.gt.1)then
693 call MPI_Reduce(newnss,g_newnss,1,
694 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
695 call MPI_Gather(newnss,1,MPI_INTEGER,
696 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
699 displ(i)=i_newnss(i-1)+displ(i-1)
701 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
702 & g_newihpb,i_newnss,displ,MPI_INTEGER,
704 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
705 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
707 if(fg_rank.eq.0) then
708 c print *,'g_newnss',g_newnss
709 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
710 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
713 newihpb(i)=g_newihpb(i)
714 newjhpb(i)=g_newjhpb(i)
722 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
727 if (idssb(i).eq.newihpb(j) .and.
728 & jdssb(i).eq.newjhpb(j)) found=.true.
732 c if (.not.found.and.fg_rank.eq.0)
733 c & write(iout,'(a15,f12.2,f8.1,2i5)')
734 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
742 if (newihpb(i).eq.idssb(j) .and.
743 & newjhpb(i).eq.jdssb(j)) found=.true.
747 c if (.not.found.and.fg_rank.eq.0)
748 c & write(iout,'(a15,f12.2,f8.1,2i5)')
749 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
763 c----------------------------------------------------------------------------
766 subroutine read_ssHist
771 include "DIMENSIONS.FREE"
772 include 'COMMON.FREE'
776 character*80 controlcard
779 call card_concat(controlcard,.true.)
781 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
788 c----------------------------------------------------------------------------
791 C-----------------------------------------------------------------------------
792 C-----------------------------------------------------------------------------
793 C-----------------------------------------------------------------------------
794 C-----------------------------------------------------------------------------
795 C-----------------------------------------------------------------------------
796 C-----------------------------------------------------------------------------
797 C-----------------------------------------------------------------------------
799 c$$$c-----------------------------------------------------------------------------
801 c$$$ subroutine ss_relax(i_in,j_in)
805 c$$$ include 'DIMENSIONS'
806 c$$$ include 'COMMON.VAR'
807 c$$$ include 'COMMON.CHAIN'
808 c$$$ include 'COMMON.IOUNITS'
809 c$$$ include 'COMMON.INTERACT'
811 c$$$c Input arguments
812 c$$$ integer i_in,j_in
814 c$$$c Local variables
815 c$$$ integer i,iretcode,nfun_sc
817 c$$$ double precision var(maxvar),e_sc,etot
824 c$$$ mask_side(i_in)=1
825 c$$$ mask_side(j_in)=1
827 c$$$c Minimize the two selected side-chains
828 c$$$ call overlap_sc(scfail) ! Better not fail!
829 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
836 c$$$c-------------------------------------------------------------
838 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
839 c$$$c Minimize side-chains only, starting from geom but without modifying
841 c$$$c If mask_r is already set, only the selected side-chains are minimized,
842 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
846 c$$$ include 'DIMENSIONS'
847 c$$$ include 'COMMON.IOUNITS'
848 c$$$ include 'COMMON.VAR'
849 c$$$ include 'COMMON.CHAIN'
850 c$$$ include 'COMMON.GEO'
851 c$$$ include 'COMMON.MINIM'
853 c$$$ common /srutu/ icall
855 c$$$c Output arguments
856 c$$$ double precision etot_sc
857 c$$$ integer iretcode,nfun
859 c$$$c External functions/subroutines
860 c$$$ external func_sc,grad_sc,fdum
862 c$$$c Local variables
864 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
866 c$$$ double precision rdum(1)
867 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
869 c$$$ integer i,nvar_restr
872 c$$$cmc start_minim=.true.
873 c$$$ call deflt(2,iv,liv,lv,v)
874 c$$$* 12 means fresh start, dont call deflt
876 c$$$* max num of fun calls
877 c$$$ if (maxfun.eq.0) maxfun=500
879 c$$$* max num of iterations
880 c$$$ if (maxmin.eq.0) maxmin=1000
882 c$$$* controls output
884 c$$$* selects output unit
886 c$$$c iv(21)=iout ! DEBUG
887 c$$$c iv(21)=8 ! DEBUG
888 c$$$* 1 means to print out result
890 c$$$c iv(22)=1 ! DEBUG
891 c$$$* 1 means to print out summary stats
893 c$$$c iv(23)=1 ! DEBUG
894 c$$$* 1 means to print initial x and d
896 c$$$c iv(24)=1 ! DEBUG
897 c$$$* min val for v(radfac) default is 0.1
899 c$$$* max val for v(radfac) default is 4.0
902 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
903 c$$$* the sumsl default is 0.1
905 c$$$* false conv if (act fnctn decrease) .lt. v(34)
906 c$$$* the sumsl default is 100*machep
907 c$$$ v(34)=v(34)/100.0D0
908 c$$$* absolute convergence
909 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
911 c$$$* relative convergence
912 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
914 c$$$* controls initial step size
916 c$$$* large vals of d correspond to small components of step
920 c$$$ do i=nphi+1,nvar
924 c$$$ call geom_to_var(nvar,x)
925 c$$$ IF (mask_r) THEN
926 c$$$ do i=1,nres ! Just in case...
930 c$$$ call x2xx(x,xx,nvar_restr)
931 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
932 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
935 c$$$c When minimizing ALL side-chains, etotal_sc is a little
936 c$$$c faster if we don't set mask_r
942 c$$$ call x2xx(x,xx,nvar_restr)
943 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
944 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
947 c$$$ call var_to_geom(nvar,x)
948 c$$$ call chainbuild_sc
955 c$$$C--------------------------------------------------------------------------
957 c$$$ subroutine chainbuild_sc
959 c$$$ include 'DIMENSIONS'
960 c$$$ include 'COMMON.VAR'
961 c$$$ include 'COMMON.INTERACT'
963 c$$$c Local variables
968 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
969 c$$$ call locate_side_chain(i)
976 c$$$C--------------------------------------------------------------------------
978 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
982 c$$$ include 'DIMENSIONS'
983 c$$$ include 'COMMON.DERIV'
984 c$$$ include 'COMMON.VAR'
985 c$$$ include 'COMMON.MINIM'
986 c$$$ include 'COMMON.IOUNITS'
988 c$$$c Input arguments
990 c$$$ double precision x(maxvar)
991 c$$$ double precision ufparm
994 c$$$c Input/Output arguments
996 c$$$ integer uiparm(1)
997 c$$$ double precision urparm(1)
999 c$$$c Output arguments
1000 c$$$ double precision f
1002 c$$$c Local variables
1003 c$$$ double precision energia(0:n_ene)
1005 c$$$c Variables used to intercept NaNs
1006 c$$$ double precision x_sum
1012 c$$$ icg=mod(nf,2)+1
1015 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
1018 c$$$ x_sum=x_sum+x(i_NAN)
1020 c$$$c Calculate the energy only if the coordinates are ok
1021 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
1022 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
1028 c$$$ call var_to_geom_restr(n,x)
1030 c$$$ call chainbuild_sc
1031 c$$$ call etotal_sc(energia(0))
1033 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
1042 c$$$c-------------------------------------------------------
1044 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
1048 c$$$ include 'DIMENSIONS'
1049 c$$$ include 'COMMON.CHAIN'
1050 c$$$ include 'COMMON.DERIV'
1051 c$$$ include 'COMMON.VAR'
1052 c$$$ include 'COMMON.INTERACT'
1053 c$$$ include 'COMMON.MINIM'
1055 c$$$c Input arguments
1057 c$$$ double precision x(maxvar)
1058 c$$$ double precision ufparm
1059 c$$$ external ufparm
1061 c$$$c Input/Output arguments
1063 c$$$ integer uiparm(1)
1064 c$$$ double precision urparm(1)
1066 c$$$c Output arguments
1067 c$$$ double precision g(maxvar)
1069 c$$$c Local variables
1070 c$$$ double precision f,gphii,gthetai,galphai,gomegai
1071 c$$$ integer ig,ind,i,j,k,igall,ij
1074 c$$$ icg=mod(nf,2)+1
1075 c$$$ if (nf-nfl+1) 20,30,40
1076 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1077 c$$$c write (iout,*) 'grad 20'
1078 c$$$ if (nf.eq.0) return
1080 c$$$ 30 call var_to_geom_restr(n,x)
1081 c$$$ call chainbuild_sc
1083 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1085 c$$$ 40 call cartder
1087 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1093 c$$$ IF (mask_phi(i+2).eq.1) THEN
1095 c$$$ do j=i+1,nres-1
1098 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1099 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1105 c$$$ ind=ind+nres-1-i
1112 c$$$ IF (mask_theta(i+2).eq.1) THEN
1115 c$$$ do j=i+1,nres-1
1118 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1119 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1124 c$$$ ind=ind+nres-1-i
1129 c$$$ if (itype(i).ne.10) then
1130 c$$$ IF (mask_side(i).eq.1) THEN
1134 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1143 c$$$ if (itype(i).ne.10) then
1144 c$$$ IF (mask_side(i).eq.1) THEN
1148 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1156 c$$$C Add the components corresponding to local energy terms.
1163 c$$$ if (mask_phi(i).eq.1) then
1165 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1171 c$$$ if (mask_theta(i).eq.1) then
1173 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1179 c$$$ if (itype(i).ne.10) then
1181 c$$$ if (mask_side(i).eq.1) then
1183 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1190 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1196 c$$$C-----------------------------------------------------------------------------
1198 c$$$ subroutine etotal_sc(energy_sc)
1202 c$$$ include 'DIMENSIONS'
1203 c$$$ include 'COMMON.VAR'
1204 c$$$ include 'COMMON.INTERACT'
1205 c$$$ include 'COMMON.DERIV'
1206 c$$$ include 'COMMON.FFIELD'
1208 c$$$c Output arguments
1209 c$$$ double precision energy_sc(0:n_ene)
1211 c$$$c Local variables
1212 c$$$ double precision evdw,escloc
1217 c$$$ energy_sc(i)=0.0D0
1220 c$$$ if (mask_r) then
1221 c$$$ call egb_sc(evdw)
1222 c$$$ call esc_sc(escloc)
1225 c$$$ call esc(escloc)
1228 c$$$ if (evdw.eq.1.0D20) then
1229 c$$$ energy_sc(0)=evdw
1231 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1233 c$$$ energy_sc(1)=evdw
1234 c$$$ energy_sc(12)=escloc
1237 c$$$C Sum up the components of the Cartesian gradient.
1241 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1248 c$$$C-----------------------------------------------------------------------------
1250 c$$$ subroutine egb_sc(evdw)
1252 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1253 c$$$C assuming the Gay-Berne potential of interaction.
1255 c$$$ implicit real*8 (a-h,o-z)
1256 c$$$ include 'DIMENSIONS'
1257 c$$$ include 'COMMON.GEO'
1258 c$$$ include 'COMMON.VAR'
1259 c$$$ include 'COMMON.LOCAL'
1260 c$$$ include 'COMMON.CHAIN'
1261 c$$$ include 'COMMON.DERIV'
1262 c$$$ include 'COMMON.NAMES'
1263 c$$$ include 'COMMON.INTERACT'
1264 c$$$ include 'COMMON.IOUNITS'
1265 c$$$ include 'COMMON.CALC'
1266 c$$$ include 'COMMON.CONTROL'
1269 c$$$ energy_dec=.false.
1270 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1273 c$$$c if (icall.eq.0) lprn=.false.
1275 c$$$ do i=iatsc_s,iatsc_e
1277 c$$$ itypi1=itype(i+1)
1281 c$$$ dxi=dc_norm(1,nres+i)
1282 c$$$ dyi=dc_norm(2,nres+i)
1283 c$$$ dzi=dc_norm(3,nres+i)
1284 c$$$c dsci_inv=dsc_inv(itypi)
1285 c$$$ dsci_inv=vbld_inv(i+nres)
1286 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1287 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1289 c$$$C Calculate SC interaction energy.
1291 c$$$ do iint=1,nint_gr(i)
1292 c$$$ do j=istart(i,iint),iend(i,iint)
1293 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1296 c$$$c dscj_inv=dsc_inv(itypj)
1297 c$$$ dscj_inv=vbld_inv(j+nres)
1298 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1299 c$$$c & 1.0d0/vbld(j+nres)
1300 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1301 c$$$ sig0ij=sigma(itypi,itypj)
1302 c$$$ chi1=chi(itypi,itypj)
1303 c$$$ chi2=chi(itypj,itypi)
1304 c$$$ chi12=chi1*chi2
1305 c$$$ chip1=chip(itypi)
1306 c$$$ chip2=chip(itypj)
1307 c$$$ chip12=chip1*chip2
1308 c$$$ alf1=alp(itypi)
1309 c$$$ alf2=alp(itypj)
1310 c$$$ alf12=0.5D0*(alf1+alf2)
1311 c$$$C For diagnostics only!!!
1321 c$$$ xj=c(1,nres+j)-xi
1322 c$$$ yj=c(2,nres+j)-yi
1323 c$$$ zj=c(3,nres+j)-zi
1324 c$$$ dxj=dc_norm(1,nres+j)
1325 c$$$ dyj=dc_norm(2,nres+j)
1326 c$$$ dzj=dc_norm(3,nres+j)
1327 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1328 c$$$c write (iout,*) "j",j," dc_norm",
1329 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1330 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1331 c$$$ rij=dsqrt(rrij)
1332 c$$$C Calculate angle-dependent terms of energy and contributions to their
1334 c$$$ call sc_angular
1335 c$$$ sigsq=1.0D0/sigsq
1336 c$$$ sig=sig0ij*dsqrt(sigsq)
1337 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1338 c$$$c for diagnostics; uncomment
1339 c$$$c rij_shift=1.2*sig0ij
1340 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1341 c$$$ if (rij_shift.le.0.0D0) then
1343 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1344 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1345 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1348 c$$$ sigder=-sig*sigsq
1349 c$$$c---------------------------------------------------------------
1350 c$$$ rij_shift=1.0D0/rij_shift
1351 c$$$ fac=rij_shift**expon
1352 c$$$ e1=fac*fac*aa(itypi,itypj)
1353 c$$$ e2=fac*bb(itypi,itypj)
1354 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1355 c$$$ eps2der=evdwij*eps3rt
1356 c$$$ eps3der=evdwij*eps2rt
1357 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1358 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1359 c$$$ evdwij=evdwij*eps2rt*eps3rt
1360 c$$$ evdw=evdw+evdwij
1362 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1363 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1364 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1365 c$$$ & restyp(itypi),i,restyp(itypj),j,
1366 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1367 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1368 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1372 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1373 c$$$ & 'evdw',i,j,evdwij
1375 c$$$C Calculate gradient components.
1376 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1377 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1378 c$$$ sigder=fac*sigder
1381 c$$$C Calculate the radial part of the gradient
1385 c$$$C Calculate angular part of the gradient.
1391 c$$$ energy_dec=.false.
1395 c$$$c-----------------------------------------------------------------------------
1397 c$$$ subroutine esc_sc(escloc)
1398 c$$$C Calculate the local energy of a side chain and its derivatives in the
1399 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1400 c$$$C ALPHA and OMEGA.
1401 c$$$ implicit real*8 (a-h,o-z)
1402 c$$$ include 'DIMENSIONS'
1403 c$$$ include 'COMMON.GEO'
1404 c$$$ include 'COMMON.LOCAL'
1405 c$$$ include 'COMMON.VAR'
1406 c$$$ include 'COMMON.INTERACT'
1407 c$$$ include 'COMMON.DERIV'
1408 c$$$ include 'COMMON.CHAIN'
1409 c$$$ include 'COMMON.IOUNITS'
1410 c$$$ include 'COMMON.NAMES'
1411 c$$$ include 'COMMON.FFIELD'
1412 c$$$ include 'COMMON.CONTROL'
1413 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1414 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1415 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1416 c$$$ delta=0.02d0*pi
1418 c$$$c write (iout,'(a)') 'ESC'
1419 c$$$ do i=loc_start,loc_end
1420 c$$$ IF (mask_side(i).eq.1) THEN
1422 c$$$ if (it.eq.10) goto 1
1423 c$$$ nlobit=nlob(it)
1424 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1425 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1426 c$$$ theti=theta(i+1)-pipol
1427 c$$$ x(1)=dtan(theti)
1431 c$$$ if (x(2).gt.pi-delta) then
1433 c$$$ xtemp(2)=pi-delta
1435 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1437 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1438 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1439 c$$$ & escloci,dersc(2))
1440 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1441 c$$$ & ddersc0(1),dersc(1))
1442 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1443 c$$$ & ddersc0(3),dersc(3))
1444 c$$$ xtemp(2)=pi-delta
1445 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1447 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1448 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1449 c$$$ & dersc0(2),esclocbi,dersc02)
1450 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1451 c$$$ & dersc12,dersc01)
1452 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1453 c$$$ dersc0(1)=dersc01
1454 c$$$ dersc0(2)=dersc02
1455 c$$$ dersc0(3)=0.0d0
1457 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1459 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1460 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1461 c$$$c & esclocbi,ss,ssd
1462 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1463 c$$$c escloci=esclocbi
1464 c$$$c write (iout,*) escloci
1465 c$$$ else if (x(2).lt.delta) then
1469 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1471 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1472 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1473 c$$$ & escloci,dersc(2))
1474 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1475 c$$$ & ddersc0(1),dersc(1))
1476 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1477 c$$$ & ddersc0(3),dersc(3))
1479 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1481 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1482 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1483 c$$$ & dersc0(2),esclocbi,dersc02)
1484 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1485 c$$$ & dersc12,dersc01)
1486 c$$$ dersc0(1)=dersc01
1487 c$$$ dersc0(2)=dersc02
1488 c$$$ dersc0(3)=0.0d0
1489 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1491 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1493 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1494 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1495 c$$$c & esclocbi,ss,ssd
1496 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1497 c$$$c write (iout,*) escloci
1499 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1502 c$$$ escloc=escloc+escloci
1503 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1504 c$$$ & 'escloc',i,escloci
1505 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1507 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1508 c$$$ & wscloc*dersc(1)
1509 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1510 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1517 c$$$C-----------------------------------------------------------------------------
1519 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1521 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1522 c$$$C assuming the Gay-Berne potential of interaction.
1524 c$$$ implicit real*8 (a-h,o-z)
1525 c$$$ include 'DIMENSIONS'
1526 c$$$ include 'COMMON.GEO'
1527 c$$$ include 'COMMON.VAR'
1528 c$$$ include 'COMMON.LOCAL'
1529 c$$$ include 'COMMON.CHAIN'
1530 c$$$ include 'COMMON.DERIV'
1531 c$$$ include 'COMMON.NAMES'
1532 c$$$ include 'COMMON.INTERACT'
1533 c$$$ include 'COMMON.IOUNITS'
1534 c$$$ include 'COMMON.CALC'
1535 c$$$ include 'COMMON.CONTROL'
1538 c$$$ energy_dec=.false.
1539 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1543 c$$$c$$$ do i=iatsc_s,iatsc_e
1546 c$$$ itypi1=itype(i+1)
1550 c$$$ dxi=dc_norm(1,nres+i)
1551 c$$$ dyi=dc_norm(2,nres+i)
1552 c$$$ dzi=dc_norm(3,nres+i)
1553 c$$$c dsci_inv=dsc_inv(itypi)
1554 c$$$ dsci_inv=vbld_inv(i+nres)
1555 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1556 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1558 c$$$C Calculate SC interaction energy.
1560 c$$$c$$$ do iint=1,nint_gr(i)
1561 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1565 c$$$c dscj_inv=dsc_inv(itypj)
1566 c$$$ dscj_inv=vbld_inv(j+nres)
1567 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1568 c$$$c & 1.0d0/vbld(j+nres)
1569 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1570 c$$$ sig0ij=sigma(itypi,itypj)
1571 c$$$ chi1=chi(itypi,itypj)
1572 c$$$ chi2=chi(itypj,itypi)
1573 c$$$ chi12=chi1*chi2
1574 c$$$ chip1=chip(itypi)
1575 c$$$ chip2=chip(itypj)
1576 c$$$ chip12=chip1*chip2
1577 c$$$ alf1=alp(itypi)
1578 c$$$ alf2=alp(itypj)
1579 c$$$ alf12=0.5D0*(alf1+alf2)
1580 c$$$C For diagnostics only!!!
1590 c$$$ xj=c(1,nres+j)-xi
1591 c$$$ yj=c(2,nres+j)-yi
1592 c$$$ zj=c(3,nres+j)-zi
1593 c$$$ dxj=dc_norm(1,nres+j)
1594 c$$$ dyj=dc_norm(2,nres+j)
1595 c$$$ dzj=dc_norm(3,nres+j)
1596 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1597 c$$$c write (iout,*) "j",j," dc_norm",
1598 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1599 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1600 c$$$ rij=dsqrt(rrij)
1601 c$$$C Calculate angle-dependent terms of energy and contributions to their
1603 c$$$ call sc_angular
1604 c$$$ sigsq=1.0D0/sigsq
1605 c$$$ sig=sig0ij*dsqrt(sigsq)
1606 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1607 c$$$c for diagnostics; uncomment
1608 c$$$c rij_shift=1.2*sig0ij
1609 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1610 c$$$ if (rij_shift.le.0.0D0) then
1612 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1613 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1614 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1617 c$$$ sigder=-sig*sigsq
1618 c$$$c---------------------------------------------------------------
1619 c$$$ rij_shift=1.0D0/rij_shift
1620 c$$$ fac=rij_shift**expon
1621 c$$$ e1=fac*fac*aa(itypi,itypj)
1622 c$$$ e2=fac*bb(itypi,itypj)
1623 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1624 c$$$ eps2der=evdwij*eps3rt
1625 c$$$ eps3der=evdwij*eps2rt
1626 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1627 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1628 c$$$ evdwij=evdwij*eps2rt*eps3rt
1629 c$$$ evdw=evdw+evdwij
1631 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1632 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1633 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1634 c$$$ & restyp(itypi),i,restyp(itypj),j,
1635 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1636 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1637 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1641 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1642 c$$$ & 'evdw',i,j,evdwij
1644 c$$$C Calculate gradient components.
1645 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1646 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1647 c$$$ sigder=fac*sigder
1650 c$$$C Calculate the radial part of the gradient
1654 c$$$C Calculate angular part of the gradient.
1657 c$$$c$$$ enddo ! iint
1659 c$$$ energy_dec=.false.
1663 c$$$C-----------------------------------------------------------------------------
1665 c$$$ subroutine perturb_side_chain(i,angle)
1669 c$$$ include 'DIMENSIONS'
1670 c$$$ include 'COMMON.CHAIN'
1671 c$$$ include 'COMMON.GEO'
1672 c$$$ include 'COMMON.VAR'
1673 c$$$ include 'COMMON.LOCAL'
1674 c$$$ include 'COMMON.IOUNITS'
1676 c$$$c External functions
1677 c$$$ external ran_number
1678 c$$$ double precision ran_number
1680 c$$$c Input arguments
1682 c$$$ double precision angle ! In degrees
1684 c$$$c Local variables
1686 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1690 c$$$ rad_ang=angle*deg2rad
1693 c$$$ do while (length.lt.0.01)
1694 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1695 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1696 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1697 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1698 c$$$ + rand_v(3)*rand_v(3)
1699 c$$$ length=sqrt(length)
1700 c$$$ rand_v(1)=rand_v(1)/length
1701 c$$$ rand_v(2)=rand_v(2)/length
1702 c$$$ rand_v(3)=rand_v(3)/length
1703 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1704 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1705 c$$$ length=1.0D0-cost*cost
1706 c$$$ if (length.lt.0.0D0) length=0.0D0
1707 c$$$ length=sqrt(length)
1708 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1709 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1710 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1712 c$$$ rand_v(1)=rand_v(1)/length
1713 c$$$ rand_v(2)=rand_v(2)/length
1714 c$$$ rand_v(3)=rand_v(3)/length
1716 c$$$ cost=dcos(rad_ang)
1717 c$$$ sint=dsin(rad_ang)
1718 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1719 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1720 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1721 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1722 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1723 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1724 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1725 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1726 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1728 c$$$ call chainbuild_cart
1733 c$$$c----------------------------------------------------------------------------
1735 c$$$ subroutine ss_relax3(i_in,j_in)
1739 c$$$ include 'DIMENSIONS'
1740 c$$$ include 'COMMON.VAR'
1741 c$$$ include 'COMMON.CHAIN'
1742 c$$$ include 'COMMON.IOUNITS'
1743 c$$$ include 'COMMON.INTERACT'
1745 c$$$c External functions
1746 c$$$ external ran_number
1747 c$$$ double precision ran_number
1749 c$$$c Input arguments
1750 c$$$ integer i_in,j_in
1752 c$$$c Local variables
1753 c$$$ double precision energy_sc(0:n_ene),etot
1754 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1755 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1756 c$$$ integer n,i_pert,i
1757 c$$$ logical notdone
1766 c$$$ mask_side(i_in)=1
1767 c$$$ mask_side(j_in)=1
1769 c$$$ call etotal_sc(energy_sc)
1770 c$$$ etot=energy_sc(0)
1771 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1772 c$$$c + energy_sc(1),energy_sc(12)
1776 c$$$ do while (notdone)
1777 c$$$ if (mod(n,2).eq.0) then
1785 c$$$ org_dc(i)=dc(i,i_pert+nres)
1786 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1787 c$$$ org_c(i)=c(i,i_pert+nres)
1789 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1790 c$$$ call perturb_side_chain(i_pert,ang_pert)
1791 c$$$ call etotal_sc(energy_sc)
1792 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1793 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1794 c$$$ if (rand_fact.lt.exp_fact) then
1795 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1796 c$$$c + energy_sc(1),energy_sc(12)
1797 c$$$ etot=energy_sc(0)
1799 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1800 c$$$c + energy_sc(1),energy_sc(12)
1802 c$$$ dc(i,i_pert+nres)=org_dc(i)
1803 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1804 c$$$ c(i,i_pert+nres)=org_c(i)
1808 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1816 c$$$c----------------------------------------------------------------------------
1818 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1820 c$$$ include 'DIMENSIONS'
1822 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1823 c$$$*********************************************************************
1824 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1825 c$$$* the calling subprogram. *
1826 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1827 c$$$* calculated in the usual pythagorean way. *
1828 c$$$* absolute convergence occurs when the function is within v(31) of *
1829 c$$$* zero. unless you know the minimum value in advance, abs convg *
1830 c$$$* is probably not useful. *
1831 c$$$* relative convergence is when the model predicts that the function *
1832 c$$$* will decrease by less than v(32)*abs(fun). *
1833 c$$$*********************************************************************
1834 c$$$ include 'COMMON.IOUNITS'
1835 c$$$ include 'COMMON.VAR'
1836 c$$$ include 'COMMON.GEO'
1837 c$$$ include 'COMMON.MINIM'
1838 c$$$ include 'COMMON.CHAIN'
1840 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1841 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1842 c$$$ + orig_ss_dist(maxres2,maxres2)
1844 c$$$ double precision etot
1845 c$$$ integer iretcode,nfun,i_in,j_in
1848 c$$$ double precision dist
1849 c$$$ external ss_func,fdum
1850 c$$$ double precision ss_func,fdum
1852 c$$$ integer iv(liv),uiparm(2)
1853 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1857 c$$$ call deflt(2,iv,liv,lv,v)
1858 c$$$* 12 means fresh start, dont call deflt
1860 c$$$* max num of fun calls
1861 c$$$ if (maxfun.eq.0) maxfun=500
1863 c$$$* max num of iterations
1864 c$$$ if (maxmin.eq.0) maxmin=1000
1866 c$$$* controls output
1868 c$$$* selects output unit
1871 c$$$* 1 means to print out result
1873 c$$$* 1 means to print out summary stats
1875 c$$$* 1 means to print initial x and d
1877 c$$$* min val for v(radfac) default is 0.1
1879 c$$$* max val for v(radfac) default is 4.0
1882 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1883 c$$$* the sumsl default is 0.1
1885 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1886 c$$$* the sumsl default is 100*machep
1887 c$$$ v(34)=v(34)/100.0D0
1888 c$$$* absolute convergence
1889 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1892 c$$$* relative convergence
1893 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1896 c$$$* controls initial step size
1898 c$$$* large vals of d correspond to small components of step
1905 c$$$ orig_ss_dc(j,i)=dc(j,i)
1908 c$$$ call geom_to_var(nvar,orig_ss_var)
1912 c$$$ orig_ss_dist(j,i)=dist(j,i)
1913 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1914 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1915 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1927 c$$$ if (ialph(i,1).gt.0) then
1930 c$$$ x(k)=dc(j,i+nres)
1937 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1940 c$$$ nfun=iv(6)+iv(30)
1950 c$$$ if (ialph(i,1).gt.0) then
1953 c$$$ dc(j,i+nres)=x(k)
1957 c$$$ call chainbuild_cart
1962 c$$$C-----------------------------------------------------------------------------
1964 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1966 c$$$ include 'DIMENSIONS'
1967 c$$$ include 'COMMON.DERIV'
1968 c$$$ include 'COMMON.IOUNITS'
1969 c$$$ include 'COMMON.VAR'
1970 c$$$ include 'COMMON.CHAIN'
1971 c$$$ include 'COMMON.INTERACT'
1972 c$$$ include 'COMMON.SBRIDGE'
1974 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1975 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1976 c$$$ + orig_ss_dist(maxres2,maxres2)
1979 c$$$ double precision x(maxres6)
1981 c$$$ double precision f
1982 c$$$ integer uiparm(2)
1983 c$$$ real*8 urparm(1)
1984 c$$$ external ufparm
1985 c$$$ double precision ufparm
1988 c$$$ double precision dist
1990 c$$$ integer i,j,k,ss_i,ss_j
1991 c$$$ double precision tempf,var(maxvar)
2006 c$$$ if (ialph(i,1).gt.0) then
2009 c$$$ dc(j,i+nres)=x(k)
2013 c$$$ call chainbuild_cart
2015 c$$$ call geom_to_var(nvar,var)
2017 c$$$c Constraints on all angles
2019 c$$$ tempf=var(i)-orig_ss_var(i)
2020 c$$$ f=f+tempf*tempf
2023 c$$$c Constraints on all distances
2025 c$$$ if (i.gt.1) then
2026 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
2027 c$$$ f=f+tempf*tempf
2030 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
2031 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
2032 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
2033 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2034 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
2035 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2036 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
2037 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2041 c$$$c Constraints for the relevant CYS-CYS
2042 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
2043 c$$$ f=f+tempf*tempf
2044 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
2046 c$$$c$$$ if (nf.ne.nfl) then
2047 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
2048 c$$$c$$$ + f,dist(5+nres,14+nres)
2056 c$$$C-----------------------------------------------------------------------------
2057 c$$$C-----------------------------------------------------------------------------
2058 subroutine triple_ssbond_ene(resi,resj,resk,eij)
2059 include 'DIMENSIONS'
2060 include 'COMMON.SBRIDGE'
2061 include 'COMMON.CHAIN'
2062 include 'COMMON.DERIV'
2063 include 'COMMON.LOCAL'
2064 include 'COMMON.INTERACT'
2065 include 'COMMON.VAR'
2066 include 'COMMON.IOUNITS'
2067 include 'COMMON.CALC'
2070 C include 'COMMON.MD'
2074 c External functions
2075 double precision h_base
2079 integer resi,resj,resk
2082 double precision eij,eij1,eij2,eij3
2086 c integer itypi,itypj,k,l
2087 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2088 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2089 double precision xik,yik,zik,xjk,yjk,zjk
2090 double precision sig0ij,ljd,sig,fac,e1,e2
2091 double precision dcosom1(3),dcosom2(3),ed
2092 double precision pom1,pom2
2093 double precision ljA,ljB,ljXs
2094 double precision d_ljB(1:3)
2095 double precision ssA,ssB,ssC,ssXs
2096 double precision ssxm,ljxm,ssm,ljm
2097 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2102 C write(iout,*) resi,resj,resk
2104 dxi=dc_norm(1,nres+i)
2105 dyi=dc_norm(2,nres+i)
2106 dzi=dc_norm(3,nres+i)
2107 dsci_inv=vbld_inv(i+nres)
2117 dxj=dc_norm(1,nres+j)
2118 dyj=dc_norm(2,nres+j)
2119 dzj=dc_norm(3,nres+j)
2120 dscj_inv=vbld_inv(j+nres)
2126 dxk=dc_norm(1,nres+k)
2127 dyk=dc_norm(2,nres+k)
2128 dzk=dc_norm(3,nres+k)
2129 dscj_inv=vbld_inv(k+nres)
2139 rrij=(xij*xij+yij*yij+zij*zij)
2140 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2141 rrik=(xik*xik+yik*yik+zik*zik)
2143 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2145 C there are three combination of distances for each trisulfide bonds
2146 C The first case the ith atom is the center
2147 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2148 C distance y is second distance the a,b,c,d are parameters derived for
2149 C this problem d parameter was set as a penalty currenlty set to 1.
2150 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2151 C second case jth atom is center
2152 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2153 C the third case kth atom is the center
2154 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2159 C write(iout,*)i,j,k,eij
2160 C The energy penalty calculated now time for the gradient part
2161 C derivative over rij
2162 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2163 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2168 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2169 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2172 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2173 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2175 C now derivative over rik
2176 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2177 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2182 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2183 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2186 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2187 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2189 C now derivative over rjk
2190 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2191 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2196 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2197 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2200 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2201 gvdwc(l,k)=gvdwc(l,k)+gg(l)