1 c----------------------------------------------------------------------------
2 subroutine check_energies
7 include 'DIMENSIONS.ZSCOPT'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SBRIDGE'
12 include 'COMMON.LOCAL'
16 double precision ran_number
20 integer i,j,k,l,lmax,p,pmax
21 double precision rmin,rmax
25 double precision wi,rij,tj,pj
49 ct wi=ran_number(0.0D0,pi)
50 c wi=ran_number(0.0D0,pi/6.0D0)
52 ct tj=ran_number(0.0D0,pi)
53 ct pj=ran_number(0.0D0,pi)
54 c pj=ran_number(0.0D0,pi/6.0D0)
58 ct rij=ran_number(rmin,rmax)
60 c(1,j)=d*sin(pj)*cos(tj)
61 c(2,j)=d*sin(pj)*sin(tj)
70 dc(k,nres+i)=c(k,nres+i)-c(k,i)
71 dc_norm(k,nres+i)=dc(k,nres+i)/d
72 dc(k,nres+j)=c(k,nres+j)-c(k,j)
73 dc_norm(k,nres+j)=dc(k,nres+j)/d
76 call dyn_ssbond_ene(i,j,eij)
85 C-----------------------------------------------------------------------------
87 subroutine dyn_ssbond_ene(resi,resj,eij)
92 include 'DIMENSIONS.ZSCOPT'
93 include 'COMMON.SBRIDGE'
94 include 'COMMON.CHAIN'
95 include 'COMMON.DERIV'
96 include 'COMMON.LOCAL'
97 include 'COMMON.INTERACT'
99 include 'COMMON.IOUNITS'
100 include 'COMMON.CALC'
103 C include 'COMMON.MD'
108 double precision h_base
119 c integer itypi,itypj,k,l
120 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
121 double precision sig0ij,ljd,sig,fac,e1,e2
122 double precision dcosom1(3),dcosom2(3),ed
123 double precision pom1,pom2
124 double precision ljA,ljB,ljXs
125 double precision d_ljB(1:3)
126 double precision ssA,ssB,ssC,ssXs
127 double precision ssxm,ljxm,ssm,ljm
128 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
129 double precision f1,f2,h1,h2,hd1,hd2
130 double precision omega,delta_inv,deltasq_inv,fac1,fac2
132 double precision xm,d_xm(1:3)
133 c-------END FIRST METHOD
134 c-------SECOND METHOD
135 c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
136 c-------END SECOND METHOD
139 logical checkstop,transgrad
140 common /sschecks/ checkstop,transgrad
142 integer icheck,nicheck,jcheck,njcheck
143 double precision echeck(-1:1),deps,ssx0,ljx0
144 c-------END TESTING CODE
151 dxi=dc_norm(1,nres+i)
152 dyi=dc_norm(2,nres+i)
153 dzi=dc_norm(3,nres+i)
154 dsci_inv=vbld_inv(i+nres)
159 if (xi.lt.0) xi=xi+boxxsize
161 if (yi.lt.0) yi=yi+boxysize
163 if (zi.lt.0) zi=zi+boxzsize
164 if ((zi.gt.bordlipbot)
165 &.and.(zi.lt.bordliptop)) then
166 C the energy transfer exist
167 if (zi.lt.buflipbot) then
168 C what fraction I am in
170 & ((zi-bordlipbot)/lipbufthick)
171 C lipbufthick is thickenes of lipid buffore
172 sslipi=sscalelip(fracinbuf)
173 ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
174 elseif (zi.gt.bufliptop) then
175 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
176 sslipi=sscalelip(fracinbuf)
177 ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
191 if (xj.lt.0) xj=xj+boxxsize
193 if (yj.lt.0) yj=yj+boxysize
195 if (zj.lt.0) zj=zj+boxzsize
196 if ((zj.gt.bordlipbot)
197 &.and.(zj.lt.bordliptop)) then
198 C the energy transfer exist
199 if (zj.lt.buflipbot) then
200 C what fraction I am in
202 & ((zj-bordlipbot)/lipbufthick)
203 C lipbufthick is thickenes of lipid buffore
204 sslipj=sscalelip(fracinbuf)
205 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
206 elseif (zj.gt.bufliptop) then
207 fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
208 sslipj=sscalelip(fracinbuf)
209 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
218 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
219 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
220 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
221 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
225 dxj=dc_norm(1,nres+j)
226 dyj=dc_norm(2,nres+j)
227 dzj=dc_norm(3,nres+j)
228 dscj_inv=vbld_inv(j+nres)
230 chi1=chi(itypi,itypj)
231 chi2=chi(itypj,itypi)
238 alf12=0.5D0*(alf1+alf2)
240 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
241 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
242 c The following are set in sc_angular
246 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
247 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
248 c om12=dxi*dxj+dyi*dyj+dzi*dzj
250 rij=1.0D0/rij ! Reset this so it makes sense
252 sig0ij=sigma(itypi,itypj)
253 sig=sig0ij*dsqrt(1.0D0/sigsq)
256 ljA=eps1*eps2rt**2*eps3rt**2
259 ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
264 deltat12=om2-om1+2.0d0
269 & +akth*(deltat1*deltat1+deltat2*deltat2)
270 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
271 ssxm=ssXs-0.5D0*ssB/ssA
274 c$$$c Some extra output
275 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
276 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
277 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
278 c$$$ if (ssx0.gt.0.0d0) then
279 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
283 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
284 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
285 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
287 c-------END TESTING CODE
290 c Stop and plot energy and derivative as a function of distance
292 ssm=ssC-0.25D0*ssB*ssB/ssA
293 ljm=-0.25D0*ljB*bb/aa
295 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
303 if (.not.checkstop) then
310 if (checkstop) rij=(ssxm-1.0d0)+
311 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
312 c-------END TESTING CODE
314 if (rij.gt.ljxm) then
317 fac=(1.0D0/ljd)**expon
320 eij=eps1*eps2rt*eps3rt*(e1+e2)
321 C write(iout,*) eij,'TU?1'
324 eij=eij*eps2rt*eps3rt
327 e1=e1*eps1*eps2rt**2*eps3rt**2
328 ed=-expon*(e1+eij)/ljd
330 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
331 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
332 eom12=eij*eps1_om12+eps2der*eps2rt_om12
333 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
334 else if (rij.lt.ssxm) then
337 eij=ssA*ssd*ssd+ssB*ssd+ssC
338 C write(iout,*) 'TU?2',ssc,ssd
339 ed=2*akcm*ssd+akct*deltat12
341 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
342 eom1=-2*akth*deltat1-pom1-om2*pom2
343 eom2= 2*akth*deltat2+pom1-om1*pom2
346 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
348 d_ssxm(1)=0.5D0*akct/ssA
352 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
353 d_ljxm(2)=d_ljxm(1)*sigsq_om2
354 d_ljxm(3)=d_ljxm(1)*sigsq_om12
355 d_ljxm(1)=d_ljxm(1)*sigsq_om1
357 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
360 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
364 ssm=ssC-0.25D0*ssB*ssB/ssA
365 d_ssm(1)=0.5D0*akct*ssB/ssA
366 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
367 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
369 f1=(rij-xm)/(ssxm-xm)
370 f2=(rij-ssxm)/(xm-ssxm)
374 C write(iout,*) eij,'TU?3'
375 delta_inv=1.0d0/(xm-ssxm)
376 deltasq_inv=delta_inv*delta_inv
378 fac1=deltasq_inv*fac*(xm-rij)
379 fac2=deltasq_inv*fac*(rij-ssxm)
380 ed=delta_inv*(Ht*hd2-ssm*hd1)
381 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
382 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
383 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
386 ljm=-0.25D0*ljB*bb/aa
387 d_ljm(1)=-0.5D0*bb/aa*ljB
388 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
389 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
391 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
392 f1=(rij-ljxm)/(xm-ljxm)
393 f2=(rij-xm)/(ljxm-xm)
397 C write(iout,*) 'TU?4',ssA
398 delta_inv=1.0d0/(ljxm-xm)
399 deltasq_inv=delta_inv*delta_inv
401 fac1=deltasq_inv*fac*(ljxm-rij)
402 fac2=deltasq_inv*fac*(rij-xm)
403 ed=delta_inv*(ljm*hd2-Ht*hd1)
404 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
405 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
406 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
408 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
410 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
416 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
417 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
418 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
420 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
421 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
422 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
423 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
426 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
428 c$$$ d_ljm(k)=ljm*d_ljB(k)
432 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
433 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
434 c$$$ d_ss(2)=akct*ssd
435 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
436 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
439 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
440 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
441 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
443 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
444 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
446 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
448 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
449 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
450 c$$$ h1=h_base(f1,hd1)
451 c$$$ h2=h_base(f2,hd2)
452 c$$$ eij=ss*h1+ljf*h2
453 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
454 c$$$ deltasq_inv=delta_inv*delta_inv
455 c$$$ fac=ljf*hd2-ss*hd1
456 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
457 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
458 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
459 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
460 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
461 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
462 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
464 c$$$ havebond=.false.
465 c$$$ if (ed.gt.0.0d0) havebond=.true.
466 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
469 write(iout,*) 'havebond',havebond
473 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
474 c write(iout,'(a15,f12.2,f8.1,2i5)')
475 c & "SSBOND_E_FORM",totT,t_bath,i,j
479 dyn_ssbond_ij(i,j)=eij
480 else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
481 dyn_ssbond_ij(i,j)=1.0d300
484 c write(iout,'(a15,f12.2,f8.1,2i5)')
485 c & "SSBOND_E_BREAK",totT,t_bath,i,j
492 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
493 & "CHECKSTOP",rij,eij,ed
498 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
505 c-------END TESTING CODE
508 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
509 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
512 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
515 gvdwx(k,i)=gvdwx(k,i)-gg(k)
516 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
517 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
518 gvdwx(k,j)=gvdwx(k,j)+gg(k)
519 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
520 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
524 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
529 gvdwc(l,i)=gvdwc(l,i)-gg(l)
530 gvdwc(l,j)=gvdwc(l,j)+gg(l)
536 C-----------------------------------------------------------------------------
538 double precision function h_base(x,deriv)
539 c A smooth function going 0->1 in range [0,1]
540 c It should NOT be called outside range [0,1], it will not work there.
547 double precision deriv
553 c Two parabolas put together. First derivative zero at extrema
554 c$$$ if (x.lt.0.5D0) then
555 c$$$ h_base=2.0D0*x*x
559 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
560 c$$$ deriv=4.0D0*deriv
563 c Third degree polynomial. First derivative zero at extrema
564 h_base=x*x*(3.0d0-2.0d0*x)
565 deriv=6.0d0*x*(1.0d0-x)
567 c Fifth degree polynomial. First and second derivatives zero at extrema
569 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
571 c$$$ deriv=deriv*deriv
572 c$$$ deriv=30.0d0*xsq*deriv
577 c----------------------------------------------------------------------------
579 subroutine dyn_set_nss
580 c Adjust nss and other relevant variables based on dyn_ssbond_ij
585 include 'DIMENSIONS.ZSCOPT'
589 include 'COMMON.SBRIDGE'
590 include 'COMMON.CHAIN'
591 include 'COMMON.IOUNITS'
592 C include 'COMMON.SETUP'
595 C include 'COMMON.MD'
600 double precision emin
602 integer diff,allflag(maxdim),allnss,
603 & allihpb(maxdim),alljhpb(maxdim),
604 & newnss,newihpb(maxdim),newjhpb(maxdim)
606 integer i_newnss(1024),displ(0:1024)
607 integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
612 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
621 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
625 if (allflag(i).eq.0 .and.
626 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
627 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
631 if (emin.lt.1.0d300) then
634 if (allflag(i).eq.0 .and.
635 & (allihpb(i).eq.allihpb(imin) .or.
636 & alljhpb(i).eq.allihpb(imin) .or.
637 & allihpb(i).eq.alljhpb(imin) .or.
638 & alljhpb(i).eq.alljhpb(imin))) then
645 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
649 if (allflag(i).eq.1) then
651 newihpb(newnss)=allihpb(i)
652 newjhpb(newnss)=alljhpb(i)
657 if (nfgtasks.gt.1)then
659 call MPI_Reduce(newnss,g_newnss,1,
660 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
661 call MPI_Gather(newnss,1,MPI_INTEGER,
662 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
665 displ(i)=i_newnss(i-1)+displ(i-1)
667 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
668 & g_newihpb,i_newnss,displ,MPI_INTEGER,
670 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
671 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
673 if(fg_rank.eq.0) then
674 c print *,'g_newnss',g_newnss
675 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
676 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
679 newihpb(i)=g_newihpb(i)
680 newjhpb(i)=g_newjhpb(i)
688 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
693 if (idssb(i).eq.newihpb(j) .and.
694 & jdssb(i).eq.newjhpb(j)) found=.true.
698 c if (.not.found.and.fg_rank.eq.0)
699 c & write(iout,'(a15,f12.2,f8.1,2i5)')
700 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
708 if (newihpb(i).eq.idssb(j) .and.
709 & newjhpb(i).eq.jdssb(j)) found=.true.
713 c if (.not.found.and.fg_rank.eq.0)
714 c & write(iout,'(a15,f12.2,f8.1,2i5)')
715 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
729 c----------------------------------------------------------------------------
733 subroutine read_ssHist
738 include "DIMENSIONS.FREE"
739 include 'COMMON.FREE'
743 character*80 controlcard
746 call card_concat(controlcard,.true.)
748 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
755 c----------------------------------------------------------------------------
758 C-----------------------------------------------------------------------------
759 C-----------------------------------------------------------------------------
760 C-----------------------------------------------------------------------------
761 C-----------------------------------------------------------------------------
762 C-----------------------------------------------------------------------------
763 C-----------------------------------------------------------------------------
764 C-----------------------------------------------------------------------------
766 c$$$c-----------------------------------------------------------------------------
768 c$$$ subroutine ss_relax(i_in,j_in)
772 c$$$ include 'DIMENSIONS'
773 c$$$ include 'COMMON.VAR'
774 c$$$ include 'COMMON.CHAIN'
775 c$$$ include 'COMMON.IOUNITS'
776 c$$$ include 'COMMON.INTERACT'
778 c$$$c Input arguments
779 c$$$ integer i_in,j_in
781 c$$$c Local variables
782 c$$$ integer i,iretcode,nfun_sc
784 c$$$ double precision var(maxvar),e_sc,etot
791 c$$$ mask_side(i_in)=1
792 c$$$ mask_side(j_in)=1
794 c$$$c Minimize the two selected side-chains
795 c$$$ call overlap_sc(scfail) ! Better not fail!
796 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
803 c$$$c-------------------------------------------------------------
805 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
806 c$$$c Minimize side-chains only, starting from geom but without modifying
808 c$$$c If mask_r is already set, only the selected side-chains are minimized,
809 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
813 c$$$ include 'DIMENSIONS'
814 c$$$ include 'COMMON.IOUNITS'
815 c$$$ include 'COMMON.VAR'
816 c$$$ include 'COMMON.CHAIN'
817 c$$$ include 'COMMON.GEO'
818 c$$$ include 'COMMON.MINIM'
820 c$$$ common /srutu/ icall
822 c$$$c Output arguments
823 c$$$ double precision etot_sc
824 c$$$ integer iretcode,nfun
826 c$$$c External functions/subroutines
827 c$$$ external func_sc,grad_sc,fdum
829 c$$$c Local variables
831 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
833 c$$$ double precision rdum(1)
834 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
836 c$$$ integer i,nvar_restr
839 c$$$cmc start_minim=.true.
840 c$$$ call deflt(2,iv,liv,lv,v)
841 c$$$* 12 means fresh start, dont call deflt
843 c$$$* max num of fun calls
844 c$$$ if (maxfun.eq.0) maxfun=500
846 c$$$* max num of iterations
847 c$$$ if (maxmin.eq.0) maxmin=1000
849 c$$$* controls output
851 c$$$* selects output unit
853 c$$$c iv(21)=iout ! DEBUG
854 c$$$c iv(21)=8 ! DEBUG
855 c$$$* 1 means to print out result
857 c$$$c iv(22)=1 ! DEBUG
858 c$$$* 1 means to print out summary stats
860 c$$$c iv(23)=1 ! DEBUG
861 c$$$* 1 means to print initial x and d
863 c$$$c iv(24)=1 ! DEBUG
864 c$$$* min val for v(radfac) default is 0.1
866 c$$$* max val for v(radfac) default is 4.0
869 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
870 c$$$* the sumsl default is 0.1
872 c$$$* false conv if (act fnctn decrease) .lt. v(34)
873 c$$$* the sumsl default is 100*machep
874 c$$$ v(34)=v(34)/100.0D0
875 c$$$* absolute convergence
876 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
878 c$$$* relative convergence
879 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
881 c$$$* controls initial step size
883 c$$$* large vals of d correspond to small components of step
887 c$$$ do i=nphi+1,nvar
891 c$$$ call geom_to_var(nvar,x)
892 c$$$ IF (mask_r) THEN
893 c$$$ do i=1,nres ! Just in case...
897 c$$$ call x2xx(x,xx,nvar_restr)
898 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
899 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
902 c$$$c When minimizing ALL side-chains, etotal_sc is a little
903 c$$$c faster if we don't set mask_r
909 c$$$ call x2xx(x,xx,nvar_restr)
910 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
911 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
914 c$$$ call var_to_geom(nvar,x)
915 c$$$ call chainbuild_sc
922 c$$$C--------------------------------------------------------------------------
924 c$$$ subroutine chainbuild_sc
926 c$$$ include 'DIMENSIONS'
927 c$$$ include 'COMMON.VAR'
928 c$$$ include 'COMMON.INTERACT'
930 c$$$c Local variables
935 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
936 c$$$ call locate_side_chain(i)
943 c$$$C--------------------------------------------------------------------------
945 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
949 c$$$ include 'DIMENSIONS'
950 c$$$ include 'COMMON.DERIV'
951 c$$$ include 'COMMON.VAR'
952 c$$$ include 'COMMON.MINIM'
953 c$$$ include 'COMMON.IOUNITS'
955 c$$$c Input arguments
957 c$$$ double precision x(maxvar)
958 c$$$ double precision ufparm
961 c$$$c Input/Output arguments
963 c$$$ integer uiparm(1)
964 c$$$ double precision urparm(1)
966 c$$$c Output arguments
967 c$$$ double precision f
969 c$$$c Local variables
970 c$$$ double precision energia(0:n_ene)
972 c$$$c Variables used to intercept NaNs
973 c$$$ double precision x_sum
982 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
985 c$$$ x_sum=x_sum+x(i_NAN)
987 c$$$c Calculate the energy only if the coordinates are ok
988 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
989 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
995 c$$$ call var_to_geom_restr(n,x)
997 c$$$ call chainbuild_sc
998 c$$$ call etotal_sc(energia(0))
1000 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
1009 c$$$c-------------------------------------------------------
1011 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
1015 c$$$ include 'DIMENSIONS'
1016 c$$$ include 'COMMON.CHAIN'
1017 c$$$ include 'COMMON.DERIV'
1018 c$$$ include 'COMMON.VAR'
1019 c$$$ include 'COMMON.INTERACT'
1020 c$$$ include 'COMMON.MINIM'
1022 c$$$c Input arguments
1024 c$$$ double precision x(maxvar)
1025 c$$$ double precision ufparm
1026 c$$$ external ufparm
1028 c$$$c Input/Output arguments
1030 c$$$ integer uiparm(1)
1031 c$$$ double precision urparm(1)
1033 c$$$c Output arguments
1034 c$$$ double precision g(maxvar)
1036 c$$$c Local variables
1037 c$$$ double precision f,gphii,gthetai,galphai,gomegai
1038 c$$$ integer ig,ind,i,j,k,igall,ij
1041 c$$$ icg=mod(nf,2)+1
1042 c$$$ if (nf-nfl+1) 20,30,40
1043 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1044 c$$$c write (iout,*) 'grad 20'
1045 c$$$ if (nf.eq.0) return
1047 c$$$ 30 call var_to_geom_restr(n,x)
1048 c$$$ call chainbuild_sc
1050 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1052 c$$$ 40 call cartder
1054 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1060 c$$$ IF (mask_phi(i+2).eq.1) THEN
1062 c$$$ do j=i+1,nres-1
1065 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1066 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1072 c$$$ ind=ind+nres-1-i
1079 c$$$ IF (mask_theta(i+2).eq.1) THEN
1082 c$$$ do j=i+1,nres-1
1085 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1086 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1091 c$$$ ind=ind+nres-1-i
1096 c$$$ if (itype(i).ne.10) then
1097 c$$$ IF (mask_side(i).eq.1) THEN
1101 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1110 c$$$ if (itype(i).ne.10) then
1111 c$$$ IF (mask_side(i).eq.1) THEN
1115 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1123 c$$$C Add the components corresponding to local energy terms.
1130 c$$$ if (mask_phi(i).eq.1) then
1132 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1138 c$$$ if (mask_theta(i).eq.1) then
1140 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1146 c$$$ if (itype(i).ne.10) then
1148 c$$$ if (mask_side(i).eq.1) then
1150 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1157 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1163 c$$$C-----------------------------------------------------------------------------
1165 c$$$ subroutine etotal_sc(energy_sc)
1169 c$$$ include 'DIMENSIONS'
1170 c$$$ include 'COMMON.VAR'
1171 c$$$ include 'COMMON.INTERACT'
1172 c$$$ include 'COMMON.DERIV'
1173 c$$$ include 'COMMON.FFIELD'
1175 c$$$c Output arguments
1176 c$$$ double precision energy_sc(0:n_ene)
1178 c$$$c Local variables
1179 c$$$ double precision evdw,escloc
1184 c$$$ energy_sc(i)=0.0D0
1187 c$$$ if (mask_r) then
1188 c$$$ call egb_sc(evdw)
1189 c$$$ call esc_sc(escloc)
1192 c$$$ call esc(escloc)
1195 c$$$ if (evdw.eq.1.0D20) then
1196 c$$$ energy_sc(0)=evdw
1198 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1200 c$$$ energy_sc(1)=evdw
1201 c$$$ energy_sc(12)=escloc
1204 c$$$C Sum up the components of the Cartesian gradient.
1208 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1215 c$$$C-----------------------------------------------------------------------------
1217 c$$$ subroutine egb_sc(evdw)
1219 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1220 c$$$C assuming the Gay-Berne potential of interaction.
1222 c$$$ implicit real*8 (a-h,o-z)
1223 c$$$ include 'DIMENSIONS'
1224 c$$$ include 'COMMON.GEO'
1225 c$$$ include 'COMMON.VAR'
1226 c$$$ include 'COMMON.LOCAL'
1227 c$$$ include 'COMMON.CHAIN'
1228 c$$$ include 'COMMON.DERIV'
1229 c$$$ include 'COMMON.NAMES'
1230 c$$$ include 'COMMON.INTERACT'
1231 c$$$ include 'COMMON.IOUNITS'
1232 c$$$ include 'COMMON.CALC'
1233 c$$$ include 'COMMON.CONTROL'
1236 c$$$ energy_dec=.false.
1237 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1240 c$$$c if (icall.eq.0) lprn=.false.
1242 c$$$ do i=iatsc_s,iatsc_e
1244 c$$$ itypi1=itype(i+1)
1248 c$$$ dxi=dc_norm(1,nres+i)
1249 c$$$ dyi=dc_norm(2,nres+i)
1250 c$$$ dzi=dc_norm(3,nres+i)
1251 c$$$c dsci_inv=dsc_inv(itypi)
1252 c$$$ dsci_inv=vbld_inv(i+nres)
1253 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1254 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1256 c$$$C Calculate SC interaction energy.
1258 c$$$ do iint=1,nint_gr(i)
1259 c$$$ do j=istart(i,iint),iend(i,iint)
1260 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1263 c$$$c dscj_inv=dsc_inv(itypj)
1264 c$$$ dscj_inv=vbld_inv(j+nres)
1265 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1266 c$$$c & 1.0d0/vbld(j+nres)
1267 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1268 c$$$ sig0ij=sigma(itypi,itypj)
1269 c$$$ chi1=chi(itypi,itypj)
1270 c$$$ chi2=chi(itypj,itypi)
1271 c$$$ chi12=chi1*chi2
1272 c$$$ chip1=chip(itypi)
1273 c$$$ chip2=chip(itypj)
1274 c$$$ chip12=chip1*chip2
1275 c$$$ alf1=alp(itypi)
1276 c$$$ alf2=alp(itypj)
1277 c$$$ alf12=0.5D0*(alf1+alf2)
1278 c$$$C For diagnostics only!!!
1288 c$$$ xj=c(1,nres+j)-xi
1289 c$$$ yj=c(2,nres+j)-yi
1290 c$$$ zj=c(3,nres+j)-zi
1291 c$$$ dxj=dc_norm(1,nres+j)
1292 c$$$ dyj=dc_norm(2,nres+j)
1293 c$$$ dzj=dc_norm(3,nres+j)
1294 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1295 c$$$c write (iout,*) "j",j," dc_norm",
1296 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1297 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1298 c$$$ rij=dsqrt(rrij)
1299 c$$$C Calculate angle-dependent terms of energy and contributions to their
1301 c$$$ call sc_angular
1302 c$$$ sigsq=1.0D0/sigsq
1303 c$$$ sig=sig0ij*dsqrt(sigsq)
1304 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1305 c$$$c for diagnostics; uncomment
1306 c$$$c rij_shift=1.2*sig0ij
1307 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1308 c$$$ if (rij_shift.le.0.0D0) then
1310 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1311 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1312 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1315 c$$$ sigder=-sig*sigsq
1316 c$$$c---------------------------------------------------------------
1317 c$$$ rij_shift=1.0D0/rij_shift
1318 c$$$ fac=rij_shift**expon
1319 c$$$ e1=fac*fac*aa(itypi,itypj)
1320 c$$$ e2=fac*bb(itypi,itypj)
1321 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1322 c$$$ eps2der=evdwij*eps3rt
1323 c$$$ eps3der=evdwij*eps2rt
1324 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1325 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1326 c$$$ evdwij=evdwij*eps2rt*eps3rt
1327 c$$$ evdw=evdw+evdwij
1329 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1330 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1331 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1332 c$$$ & restyp(itypi),i,restyp(itypj),j,
1333 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1334 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1335 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1339 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1340 c$$$ & 'evdw',i,j,evdwij
1342 c$$$C Calculate gradient components.
1343 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1344 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1345 c$$$ sigder=fac*sigder
1348 c$$$C Calculate the radial part of the gradient
1352 c$$$C Calculate angular part of the gradient.
1358 c$$$ energy_dec=.false.
1362 c$$$c-----------------------------------------------------------------------------
1364 c$$$ subroutine esc_sc(escloc)
1365 c$$$C Calculate the local energy of a side chain and its derivatives in the
1366 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1367 c$$$C ALPHA and OMEGA.
1368 c$$$ implicit real*8 (a-h,o-z)
1369 c$$$ include 'DIMENSIONS'
1370 c$$$ include 'COMMON.GEO'
1371 c$$$ include 'COMMON.LOCAL'
1372 c$$$ include 'COMMON.VAR'
1373 c$$$ include 'COMMON.INTERACT'
1374 c$$$ include 'COMMON.DERIV'
1375 c$$$ include 'COMMON.CHAIN'
1376 c$$$ include 'COMMON.IOUNITS'
1377 c$$$ include 'COMMON.NAMES'
1378 c$$$ include 'COMMON.FFIELD'
1379 c$$$ include 'COMMON.CONTROL'
1380 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1381 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1382 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1383 c$$$ delta=0.02d0*pi
1385 c$$$c write (iout,'(a)') 'ESC'
1386 c$$$ do i=loc_start,loc_end
1387 c$$$ IF (mask_side(i).eq.1) THEN
1389 c$$$ if (it.eq.10) goto 1
1390 c$$$ nlobit=nlob(it)
1391 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1392 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1393 c$$$ theti=theta(i+1)-pipol
1394 c$$$ x(1)=dtan(theti)
1398 c$$$ if (x(2).gt.pi-delta) then
1400 c$$$ xtemp(2)=pi-delta
1402 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1404 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1405 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1406 c$$$ & escloci,dersc(2))
1407 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1408 c$$$ & ddersc0(1),dersc(1))
1409 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1410 c$$$ & ddersc0(3),dersc(3))
1411 c$$$ xtemp(2)=pi-delta
1412 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1414 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1415 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1416 c$$$ & dersc0(2),esclocbi,dersc02)
1417 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1418 c$$$ & dersc12,dersc01)
1419 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1420 c$$$ dersc0(1)=dersc01
1421 c$$$ dersc0(2)=dersc02
1422 c$$$ dersc0(3)=0.0d0
1424 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1426 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1427 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1428 c$$$c & esclocbi,ss,ssd
1429 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1430 c$$$c escloci=esclocbi
1431 c$$$c write (iout,*) escloci
1432 c$$$ else if (x(2).lt.delta) then
1436 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1438 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1439 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1440 c$$$ & escloci,dersc(2))
1441 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1442 c$$$ & ddersc0(1),dersc(1))
1443 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1444 c$$$ & ddersc0(3),dersc(3))
1446 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1448 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1449 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1450 c$$$ & dersc0(2),esclocbi,dersc02)
1451 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1452 c$$$ & dersc12,dersc01)
1453 c$$$ dersc0(1)=dersc01
1454 c$$$ dersc0(2)=dersc02
1455 c$$$ dersc0(3)=0.0d0
1456 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1458 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1460 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1461 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1462 c$$$c & esclocbi,ss,ssd
1463 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1464 c$$$c write (iout,*) escloci
1466 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1469 c$$$ escloc=escloc+escloci
1470 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1471 c$$$ & 'escloc',i,escloci
1472 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1474 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1475 c$$$ & wscloc*dersc(1)
1476 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1477 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1484 c$$$C-----------------------------------------------------------------------------
1486 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1488 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1489 c$$$C assuming the Gay-Berne potential of interaction.
1491 c$$$ implicit real*8 (a-h,o-z)
1492 c$$$ include 'DIMENSIONS'
1493 c$$$ include 'COMMON.GEO'
1494 c$$$ include 'COMMON.VAR'
1495 c$$$ include 'COMMON.LOCAL'
1496 c$$$ include 'COMMON.CHAIN'
1497 c$$$ include 'COMMON.DERIV'
1498 c$$$ include 'COMMON.NAMES'
1499 c$$$ include 'COMMON.INTERACT'
1500 c$$$ include 'COMMON.IOUNITS'
1501 c$$$ include 'COMMON.CALC'
1502 c$$$ include 'COMMON.CONTROL'
1505 c$$$ energy_dec=.false.
1506 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1510 c$$$c$$$ do i=iatsc_s,iatsc_e
1513 c$$$ itypi1=itype(i+1)
1517 c$$$ dxi=dc_norm(1,nres+i)
1518 c$$$ dyi=dc_norm(2,nres+i)
1519 c$$$ dzi=dc_norm(3,nres+i)
1520 c$$$c dsci_inv=dsc_inv(itypi)
1521 c$$$ dsci_inv=vbld_inv(i+nres)
1522 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1523 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1525 c$$$C Calculate SC interaction energy.
1527 c$$$c$$$ do iint=1,nint_gr(i)
1528 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1532 c$$$c dscj_inv=dsc_inv(itypj)
1533 c$$$ dscj_inv=vbld_inv(j+nres)
1534 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1535 c$$$c & 1.0d0/vbld(j+nres)
1536 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1537 c$$$ sig0ij=sigma(itypi,itypj)
1538 c$$$ chi1=chi(itypi,itypj)
1539 c$$$ chi2=chi(itypj,itypi)
1540 c$$$ chi12=chi1*chi2
1541 c$$$ chip1=chip(itypi)
1542 c$$$ chip2=chip(itypj)
1543 c$$$ chip12=chip1*chip2
1544 c$$$ alf1=alp(itypi)
1545 c$$$ alf2=alp(itypj)
1546 c$$$ alf12=0.5D0*(alf1+alf2)
1547 c$$$C For diagnostics only!!!
1557 c$$$ xj=c(1,nres+j)-xi
1558 c$$$ yj=c(2,nres+j)-yi
1559 c$$$ zj=c(3,nres+j)-zi
1560 c$$$ dxj=dc_norm(1,nres+j)
1561 c$$$ dyj=dc_norm(2,nres+j)
1562 c$$$ dzj=dc_norm(3,nres+j)
1563 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1564 c$$$c write (iout,*) "j",j," dc_norm",
1565 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1566 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1567 c$$$ rij=dsqrt(rrij)
1568 c$$$C Calculate angle-dependent terms of energy and contributions to their
1570 c$$$ call sc_angular
1571 c$$$ sigsq=1.0D0/sigsq
1572 c$$$ sig=sig0ij*dsqrt(sigsq)
1573 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1574 c$$$c for diagnostics; uncomment
1575 c$$$c rij_shift=1.2*sig0ij
1576 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1577 c$$$ if (rij_shift.le.0.0D0) then
1579 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1580 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1581 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1584 c$$$ sigder=-sig*sigsq
1585 c$$$c---------------------------------------------------------------
1586 c$$$ rij_shift=1.0D0/rij_shift
1587 c$$$ fac=rij_shift**expon
1588 c$$$ e1=fac*fac*aa(itypi,itypj)
1589 c$$$ e2=fac*bb(itypi,itypj)
1590 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1591 c$$$ eps2der=evdwij*eps3rt
1592 c$$$ eps3der=evdwij*eps2rt
1593 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1594 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1595 c$$$ evdwij=evdwij*eps2rt*eps3rt
1596 c$$$ evdw=evdw+evdwij
1598 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1599 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1600 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1601 c$$$ & restyp(itypi),i,restyp(itypj),j,
1602 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1603 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1604 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1608 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1609 c$$$ & 'evdw',i,j,evdwij
1611 c$$$C Calculate gradient components.
1612 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1613 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1614 c$$$ sigder=fac*sigder
1617 c$$$C Calculate the radial part of the gradient
1621 c$$$C Calculate angular part of the gradient.
1624 c$$$c$$$ enddo ! iint
1626 c$$$ energy_dec=.false.
1630 c$$$C-----------------------------------------------------------------------------
1632 c$$$ subroutine perturb_side_chain(i,angle)
1636 c$$$ include 'DIMENSIONS'
1637 c$$$ include 'COMMON.CHAIN'
1638 c$$$ include 'COMMON.GEO'
1639 c$$$ include 'COMMON.VAR'
1640 c$$$ include 'COMMON.LOCAL'
1641 c$$$ include 'COMMON.IOUNITS'
1643 c$$$c External functions
1644 c$$$ external ran_number
1645 c$$$ double precision ran_number
1647 c$$$c Input arguments
1649 c$$$ double precision angle ! In degrees
1651 c$$$c Local variables
1653 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1657 c$$$ rad_ang=angle*deg2rad
1660 c$$$ do while (length.lt.0.01)
1661 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1662 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1663 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1664 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1665 c$$$ + rand_v(3)*rand_v(3)
1666 c$$$ length=sqrt(length)
1667 c$$$ rand_v(1)=rand_v(1)/length
1668 c$$$ rand_v(2)=rand_v(2)/length
1669 c$$$ rand_v(3)=rand_v(3)/length
1670 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1671 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1672 c$$$ length=1.0D0-cost*cost
1673 c$$$ if (length.lt.0.0D0) length=0.0D0
1674 c$$$ length=sqrt(length)
1675 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1676 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1677 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1679 c$$$ rand_v(1)=rand_v(1)/length
1680 c$$$ rand_v(2)=rand_v(2)/length
1681 c$$$ rand_v(3)=rand_v(3)/length
1683 c$$$ cost=dcos(rad_ang)
1684 c$$$ sint=dsin(rad_ang)
1685 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1686 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1687 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1688 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1689 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1690 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1691 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1692 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1693 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1695 c$$$ call chainbuild_cart
1700 c$$$c----------------------------------------------------------------------------
1702 c$$$ subroutine ss_relax3(i_in,j_in)
1706 c$$$ include 'DIMENSIONS'
1707 c$$$ include 'COMMON.VAR'
1708 c$$$ include 'COMMON.CHAIN'
1709 c$$$ include 'COMMON.IOUNITS'
1710 c$$$ include 'COMMON.INTERACT'
1712 c$$$c External functions
1713 c$$$ external ran_number
1714 c$$$ double precision ran_number
1716 c$$$c Input arguments
1717 c$$$ integer i_in,j_in
1719 c$$$c Local variables
1720 c$$$ double precision energy_sc(0:n_ene),etot
1721 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1722 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1723 c$$$ integer n,i_pert,i
1724 c$$$ logical notdone
1733 c$$$ mask_side(i_in)=1
1734 c$$$ mask_side(j_in)=1
1736 c$$$ call etotal_sc(energy_sc)
1737 c$$$ etot=energy_sc(0)
1738 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1739 c$$$c + energy_sc(1),energy_sc(12)
1743 c$$$ do while (notdone)
1744 c$$$ if (mod(n,2).eq.0) then
1752 c$$$ org_dc(i)=dc(i,i_pert+nres)
1753 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1754 c$$$ org_c(i)=c(i,i_pert+nres)
1756 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1757 c$$$ call perturb_side_chain(i_pert,ang_pert)
1758 c$$$ call etotal_sc(energy_sc)
1759 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1760 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1761 c$$$ if (rand_fact.lt.exp_fact) then
1762 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1763 c$$$c + energy_sc(1),energy_sc(12)
1764 c$$$ etot=energy_sc(0)
1766 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1767 c$$$c + energy_sc(1),energy_sc(12)
1769 c$$$ dc(i,i_pert+nres)=org_dc(i)
1770 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1771 c$$$ c(i,i_pert+nres)=org_c(i)
1775 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1783 c$$$c----------------------------------------------------------------------------
1785 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1787 c$$$ include 'DIMENSIONS'
1789 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1790 c$$$*********************************************************************
1791 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1792 c$$$* the calling subprogram. *
1793 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1794 c$$$* calculated in the usual pythagorean way. *
1795 c$$$* absolute convergence occurs when the function is within v(31) of *
1796 c$$$* zero. unless you know the minimum value in advance, abs convg *
1797 c$$$* is probably not useful. *
1798 c$$$* relative convergence is when the model predicts that the function *
1799 c$$$* will decrease by less than v(32)*abs(fun). *
1800 c$$$*********************************************************************
1801 c$$$ include 'COMMON.IOUNITS'
1802 c$$$ include 'COMMON.VAR'
1803 c$$$ include 'COMMON.GEO'
1804 c$$$ include 'COMMON.MINIM'
1805 c$$$ include 'COMMON.CHAIN'
1807 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1808 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1809 c$$$ + orig_ss_dist(maxres2,maxres2)
1811 c$$$ double precision etot
1812 c$$$ integer iretcode,nfun,i_in,j_in
1815 c$$$ double precision dist
1816 c$$$ external ss_func,fdum
1817 c$$$ double precision ss_func,fdum
1819 c$$$ integer iv(liv),uiparm(2)
1820 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1824 c$$$ call deflt(2,iv,liv,lv,v)
1825 c$$$* 12 means fresh start, dont call deflt
1827 c$$$* max num of fun calls
1828 c$$$ if (maxfun.eq.0) maxfun=500
1830 c$$$* max num of iterations
1831 c$$$ if (maxmin.eq.0) maxmin=1000
1833 c$$$* controls output
1835 c$$$* selects output unit
1838 c$$$* 1 means to print out result
1840 c$$$* 1 means to print out summary stats
1842 c$$$* 1 means to print initial x and d
1844 c$$$* min val for v(radfac) default is 0.1
1846 c$$$* max val for v(radfac) default is 4.0
1849 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1850 c$$$* the sumsl default is 0.1
1852 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1853 c$$$* the sumsl default is 100*machep
1854 c$$$ v(34)=v(34)/100.0D0
1855 c$$$* absolute convergence
1856 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1859 c$$$* relative convergence
1860 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1863 c$$$* controls initial step size
1865 c$$$* large vals of d correspond to small components of step
1872 c$$$ orig_ss_dc(j,i)=dc(j,i)
1875 c$$$ call geom_to_var(nvar,orig_ss_var)
1879 c$$$ orig_ss_dist(j,i)=dist(j,i)
1880 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1881 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1882 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1894 c$$$ if (ialph(i,1).gt.0) then
1897 c$$$ x(k)=dc(j,i+nres)
1904 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1907 c$$$ nfun=iv(6)+iv(30)
1917 c$$$ if (ialph(i,1).gt.0) then
1920 c$$$ dc(j,i+nres)=x(k)
1924 c$$$ call chainbuild_cart
1929 c$$$C-----------------------------------------------------------------------------
1931 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1933 c$$$ include 'DIMENSIONS'
1934 c$$$ include 'COMMON.DERIV'
1935 c$$$ include 'COMMON.IOUNITS'
1936 c$$$ include 'COMMON.VAR'
1937 c$$$ include 'COMMON.CHAIN'
1938 c$$$ include 'COMMON.INTERACT'
1939 c$$$ include 'COMMON.SBRIDGE'
1941 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1942 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1943 c$$$ + orig_ss_dist(maxres2,maxres2)
1946 c$$$ double precision x(maxres6)
1948 c$$$ double precision f
1949 c$$$ integer uiparm(2)
1950 c$$$ real*8 urparm(1)
1951 c$$$ external ufparm
1952 c$$$ double precision ufparm
1955 c$$$ double precision dist
1957 c$$$ integer i,j,k,ss_i,ss_j
1958 c$$$ double precision tempf,var(maxvar)
1973 c$$$ if (ialph(i,1).gt.0) then
1976 c$$$ dc(j,i+nres)=x(k)
1980 c$$$ call chainbuild_cart
1982 c$$$ call geom_to_var(nvar,var)
1984 c$$$c Constraints on all angles
1986 c$$$ tempf=var(i)-orig_ss_var(i)
1987 c$$$ f=f+tempf*tempf
1990 c$$$c Constraints on all distances
1992 c$$$ if (i.gt.1) then
1993 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1994 c$$$ f=f+tempf*tempf
1997 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1998 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1999 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
2000 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2001 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
2002 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2003 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
2004 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2008 c$$$c Constraints for the relevant CYS-CYS
2009 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
2010 c$$$ f=f+tempf*tempf
2011 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
2013 c$$$c$$$ if (nf.ne.nfl) then
2014 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
2015 c$$$c$$$ + f,dist(5+nres,14+nres)
2023 c$$$C-----------------------------------------------------------------------------
2024 c$$$C-----------------------------------------------------------------------------
2025 subroutine triple_ssbond_ene(resi,resj,resk,eij)
2026 include 'DIMENSIONS'
2027 include 'COMMON.SBRIDGE'
2028 include 'COMMON.CHAIN'
2029 include 'COMMON.DERIV'
2030 include 'COMMON.LOCAL'
2031 include 'COMMON.INTERACT'
2032 include 'COMMON.VAR'
2033 include 'COMMON.IOUNITS'
2034 include 'COMMON.CALC'
2037 C include 'COMMON.MD'
2041 c External functions
2042 double precision h_base
2046 integer resi,resj,resk
2049 double precision eij,eij1,eij2,eij3
2053 c integer itypi,itypj,k,l
2054 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2055 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2056 double precision xik,yik,zik,xjk,yjk,zjk
2057 double precision sig0ij,ljd,sig,fac,e1,e2
2058 double precision dcosom1(3),dcosom2(3),ed
2059 double precision pom1,pom2
2060 double precision ljA,ljB,ljXs
2061 double precision d_ljB(1:3)
2062 double precision ssA,ssB,ssC,ssXs
2063 double precision ssxm,ljxm,ssm,ljm
2064 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2069 C write(iout,*) resi,resj,resk
2071 dxi=dc_norm(1,nres+i)
2072 dyi=dc_norm(2,nres+i)
2073 dzi=dc_norm(3,nres+i)
2074 dsci_inv=vbld_inv(i+nres)
2084 dxj=dc_norm(1,nres+j)
2085 dyj=dc_norm(2,nres+j)
2086 dzj=dc_norm(3,nres+j)
2087 dscj_inv=vbld_inv(j+nres)
2093 dxk=dc_norm(1,nres+k)
2094 dyk=dc_norm(2,nres+k)
2095 dzk=dc_norm(3,nres+k)
2096 dscj_inv=vbld_inv(k+nres)
2106 rrij=(xij*xij+yij*yij+zij*zij)
2107 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2108 rrik=(xik*xik+yik*yik+zik*zik)
2110 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2112 C there are three combination of distances for each trisulfide bonds
2113 C The first case the ith atom is the center
2114 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2115 C distance y is second distance the a,b,c,d are parameters derived for
2116 C this problem d parameter was set as a penalty currenlty set to 1.
2117 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2118 C second case jth atom is center
2119 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2120 C the third case kth atom is the center
2121 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2126 C write(iout,*)i,j,k,eij
2127 C The energy penalty calculated now time for the gradient part
2128 C derivative over rij
2129 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2130 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2135 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2136 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2139 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2140 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2142 C now derivative over rik
2143 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2144 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2149 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2150 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2153 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2154 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2156 C now derivative over rjk
2157 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2158 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2163 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2164 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2167 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2168 gvdwc(l,k)=gvdwc(l,k)+gg(l)