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)
157 if (xi.lt.0) xi=xi+boxxsize
159 if (yi.lt.0) yi=yi+boxysize
161 if (zi.lt.0) zi=zi+boxzsize
162 if ((zi.gt.bordlipbot)
163 &.and.(zi.lt.bordliptop)) then
164 C the energy transfer exist
165 if (zi.lt.buflipbot) then
166 C what fraction I am in
168 & ((zi-bordlipbot)/lipbufthick)
169 C lipbufthick is thickenes of lipid buffore
170 sslipi=sscalelip(fracinbuf)
171 ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
172 elseif (zi.gt.bufliptop) then
173 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
174 sslipi=sscalelip(fracinbuf)
175 ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
189 if (xj.lt.0) xj=xj+boxxsize
191 if (yj.lt.0) yj=yj+boxysize
193 if (zj.lt.0) zj=zj+boxzsize
194 if ((zj.gt.bordlipbot)
195 &.and.(zj.lt.bordliptop)) then
196 C the energy transfer exist
197 if (zj.lt.buflipbot) then
198 C what fraction I am in
200 & ((zj-bordlipbot)/lipbufthick)
201 C lipbufthick is thickenes of lipid buffore
202 sslipj=sscalelip(fracinbuf)
203 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
204 elseif (zj.gt.bufliptop) then
205 fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
206 sslipj=sscalelip(fracinbuf)
207 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
216 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
217 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
218 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
219 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
223 dxj=dc_norm(1,nres+j)
224 dyj=dc_norm(2,nres+j)
225 dzj=dc_norm(3,nres+j)
226 dscj_inv=vbld_inv(j+nres)
228 chi1=chi(itypi,itypj)
229 chi2=chi(itypj,itypi)
236 alf12=0.5D0*(alf1+alf2)
238 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
239 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
240 c The following are set in sc_angular
244 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
245 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
246 c om12=dxi*dxj+dyi*dyj+dzi*dzj
248 rij=1.0D0/rij ! Reset this so it makes sense
250 sig0ij=sigma(itypi,itypj)
251 sig=sig0ij*dsqrt(1.0D0/sigsq)
254 ljA=eps1*eps2rt**2*eps3rt**2
257 ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
262 deltat12=om2-om1+2.0d0
267 & +akth*(deltat1*deltat1+deltat2*deltat2)
268 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
269 ssxm=ssXs-0.5D0*ssB/ssA
272 c$$$c Some extra output
273 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
274 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
275 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
276 c$$$ if (ssx0.gt.0.0d0) then
277 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
281 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
282 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
283 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
285 c-------END TESTING CODE
288 c Stop and plot energy and derivative as a function of distance
290 ssm=ssC-0.25D0*ssB*ssB/ssA
291 ljm=-0.25D0*ljB*bb/aa
293 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
301 if (.not.checkstop) then
308 if (checkstop) rij=(ssxm-1.0d0)+
309 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
310 c-------END TESTING CODE
312 if (rij.gt.ljxm) then
315 fac=(1.0D0/ljd)**expon
318 eij=eps1*eps2rt*eps3rt*(e1+e2)
319 C write(iout,*) eij,'TU?1'
322 eij=eij*eps2rt*eps3rt
325 e1=e1*eps1*eps2rt**2*eps3rt**2
326 ed=-expon*(e1+eij)/ljd
328 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
329 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
330 eom12=eij*eps1_om12+eps2der*eps2rt_om12
331 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
332 else if (rij.lt.ssxm) then
335 eij=ssA*ssd*ssd+ssB*ssd+ssC
336 C write(iout,*) 'TU?2',ssc,ssd
337 ed=2*akcm*ssd+akct*deltat12
339 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
340 eom1=-2*akth*deltat1-pom1-om2*pom2
341 eom2= 2*akth*deltat2+pom1-om1*pom2
344 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
346 d_ssxm(1)=0.5D0*akct/ssA
350 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
351 d_ljxm(2)=d_ljxm(1)*sigsq_om2
352 d_ljxm(3)=d_ljxm(1)*sigsq_om12
353 d_ljxm(1)=d_ljxm(1)*sigsq_om1
355 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
358 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
362 ssm=ssC-0.25D0*ssB*ssB/ssA
363 d_ssm(1)=0.5D0*akct*ssB/ssA
364 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
365 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
367 f1=(rij-xm)/(ssxm-xm)
368 f2=(rij-ssxm)/(xm-ssxm)
372 C write(iout,*) eij,'TU?3'
373 delta_inv=1.0d0/(xm-ssxm)
374 deltasq_inv=delta_inv*delta_inv
376 fac1=deltasq_inv*fac*(xm-rij)
377 fac2=deltasq_inv*fac*(rij-ssxm)
378 ed=delta_inv*(Ht*hd2-ssm*hd1)
379 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
380 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
381 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
384 ljm=-0.25D0*ljB*bb/aa
385 d_ljm(1)=-0.5D0*bb/aa*ljB
386 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
387 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
389 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
390 f1=(rij-ljxm)/(xm-ljxm)
391 f2=(rij-xm)/(ljxm-xm)
395 C write(iout,*) 'TU?4',ssA
396 delta_inv=1.0d0/(ljxm-xm)
397 deltasq_inv=delta_inv*delta_inv
399 fac1=deltasq_inv*fac*(ljxm-rij)
400 fac2=deltasq_inv*fac*(rij-xm)
401 ed=delta_inv*(ljm*hd2-Ht*hd1)
402 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
403 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
404 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
406 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
408 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
414 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
415 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
416 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
418 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
419 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
420 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
421 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
424 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
426 c$$$ d_ljm(k)=ljm*d_ljB(k)
430 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
431 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
432 c$$$ d_ss(2)=akct*ssd
433 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
434 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
437 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
438 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
439 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
441 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
442 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
444 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
446 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
447 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
448 c$$$ h1=h_base(f1,hd1)
449 c$$$ h2=h_base(f2,hd2)
450 c$$$ eij=ss*h1+ljf*h2
451 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
452 c$$$ deltasq_inv=delta_inv*delta_inv
453 c$$$ fac=ljf*hd2-ss*hd1
454 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
455 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
456 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
457 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
458 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
459 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
460 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
462 c$$$ havebond=.false.
463 c$$$ if (ed.gt.0.0d0) havebond=.true.
464 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
467 write(iout,*) 'havebond',havebond
471 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
472 c write(iout,'(a15,f12.2,f8.1,2i5)')
473 c & "SSBOND_E_FORM",totT,t_bath,i,j
477 dyn_ssbond_ij(i,j)=eij
478 else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
479 dyn_ssbond_ij(i,j)=1.0d300
482 c write(iout,'(a15,f12.2,f8.1,2i5)')
483 c & "SSBOND_E_BREAK",totT,t_bath,i,j
490 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
491 & "CHECKSTOP",rij,eij,ed
496 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
503 c-------END TESTING CODE
506 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
507 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
510 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
513 gvdwx(k,i)=gvdwx(k,i)-gg(k)
514 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
515 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
516 gvdwx(k,j)=gvdwx(k,j)+gg(k)
517 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
518 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
522 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
527 gvdwc(l,i)=gvdwc(l,i)-gg(l)
528 gvdwc(l,j)=gvdwc(l,j)+gg(l)
534 C-----------------------------------------------------------------------------
536 double precision function h_base(x,deriv)
537 c A smooth function going 0->1 in range [0,1]
538 c It should NOT be called outside range [0,1], it will not work there.
545 double precision deriv
551 c Two parabolas put together. First derivative zero at extrema
552 c$$$ if (x.lt.0.5D0) then
553 c$$$ h_base=2.0D0*x*x
557 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
558 c$$$ deriv=4.0D0*deriv
561 c Third degree polynomial. First derivative zero at extrema
562 h_base=x*x*(3.0d0-2.0d0*x)
563 deriv=6.0d0*x*(1.0d0-x)
565 c Fifth degree polynomial. First and second derivatives zero at extrema
567 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
569 c$$$ deriv=deriv*deriv
570 c$$$ deriv=30.0d0*xsq*deriv
575 c----------------------------------------------------------------------------
577 subroutine dyn_set_nss
578 c Adjust nss and other relevant variables based on dyn_ssbond_ij
586 include 'COMMON.SBRIDGE'
587 include 'COMMON.CHAIN'
588 include 'COMMON.IOUNITS'
589 C include 'COMMON.SETUP'
592 C include 'COMMON.MD'
597 double precision emin
599 integer diff,allflag(maxdim),allnss,
600 & allihpb(maxdim),alljhpb(maxdim),
601 & newnss,newihpb(maxdim),newjhpb(maxdim)
603 integer i_newnss(1024),displ(0:1024)
604 integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
609 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
618 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
622 if (allflag(i).eq.0 .and.
623 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
624 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
628 if (emin.lt.1.0d300) then
631 if (allflag(i).eq.0 .and.
632 & (allihpb(i).eq.allihpb(imin) .or.
633 & alljhpb(i).eq.allihpb(imin) .or.
634 & allihpb(i).eq.alljhpb(imin) .or.
635 & alljhpb(i).eq.alljhpb(imin))) then
642 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
646 if (allflag(i).eq.1) then
648 newihpb(newnss)=allihpb(i)
649 newjhpb(newnss)=alljhpb(i)
654 if (nfgtasks.gt.1)then
656 call MPI_Reduce(newnss,g_newnss,1,
657 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
658 call MPI_Gather(newnss,1,MPI_INTEGER,
659 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
662 displ(i)=i_newnss(i-1)+displ(i-1)
664 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
665 & g_newihpb,i_newnss,displ,MPI_INTEGER,
667 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
668 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
670 if(fg_rank.eq.0) then
671 c print *,'g_newnss',g_newnss
672 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
673 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
676 newihpb(i)=g_newihpb(i)
677 newjhpb(i)=g_newjhpb(i)
685 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
690 if (idssb(i).eq.newihpb(j) .and.
691 & jdssb(i).eq.newjhpb(j)) found=.true.
695 c if (.not.found.and.fg_rank.eq.0)
696 c & write(iout,'(a15,f12.2,f8.1,2i5)')
697 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
705 if (newihpb(i).eq.idssb(j) .and.
706 & newjhpb(i).eq.jdssb(j)) found=.true.
710 c if (.not.found.and.fg_rank.eq.0)
711 c & write(iout,'(a15,f12.2,f8.1,2i5)')
712 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
726 c----------------------------------------------------------------------------
729 subroutine read_ssHist
734 include "DIMENSIONS.FREE"
735 include 'COMMON.FREE'
739 character*80 controlcard
742 call card_concat(controlcard,.true.)
744 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
751 c----------------------------------------------------------------------------
754 C-----------------------------------------------------------------------------
755 C-----------------------------------------------------------------------------
756 C-----------------------------------------------------------------------------
757 C-----------------------------------------------------------------------------
758 C-----------------------------------------------------------------------------
759 C-----------------------------------------------------------------------------
760 C-----------------------------------------------------------------------------
762 c$$$c-----------------------------------------------------------------------------
764 c$$$ subroutine ss_relax(i_in,j_in)
768 c$$$ include 'DIMENSIONS'
769 c$$$ include 'COMMON.VAR'
770 c$$$ include 'COMMON.CHAIN'
771 c$$$ include 'COMMON.IOUNITS'
772 c$$$ include 'COMMON.INTERACT'
774 c$$$c Input arguments
775 c$$$ integer i_in,j_in
777 c$$$c Local variables
778 c$$$ integer i,iretcode,nfun_sc
780 c$$$ double precision var(maxvar),e_sc,etot
787 c$$$ mask_side(i_in)=1
788 c$$$ mask_side(j_in)=1
790 c$$$c Minimize the two selected side-chains
791 c$$$ call overlap_sc(scfail) ! Better not fail!
792 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
799 c$$$c-------------------------------------------------------------
801 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
802 c$$$c Minimize side-chains only, starting from geom but without modifying
804 c$$$c If mask_r is already set, only the selected side-chains are minimized,
805 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
809 c$$$ include 'DIMENSIONS'
810 c$$$ include 'COMMON.IOUNITS'
811 c$$$ include 'COMMON.VAR'
812 c$$$ include 'COMMON.CHAIN'
813 c$$$ include 'COMMON.GEO'
814 c$$$ include 'COMMON.MINIM'
816 c$$$ common /srutu/ icall
818 c$$$c Output arguments
819 c$$$ double precision etot_sc
820 c$$$ integer iretcode,nfun
822 c$$$c External functions/subroutines
823 c$$$ external func_sc,grad_sc,fdum
825 c$$$c Local variables
827 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
829 c$$$ double precision rdum(1)
830 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
832 c$$$ integer i,nvar_restr
835 c$$$cmc start_minim=.true.
836 c$$$ call deflt(2,iv,liv,lv,v)
837 c$$$* 12 means fresh start, dont call deflt
839 c$$$* max num of fun calls
840 c$$$ if (maxfun.eq.0) maxfun=500
842 c$$$* max num of iterations
843 c$$$ if (maxmin.eq.0) maxmin=1000
845 c$$$* controls output
847 c$$$* selects output unit
849 c$$$c iv(21)=iout ! DEBUG
850 c$$$c iv(21)=8 ! DEBUG
851 c$$$* 1 means to print out result
853 c$$$c iv(22)=1 ! DEBUG
854 c$$$* 1 means to print out summary stats
856 c$$$c iv(23)=1 ! DEBUG
857 c$$$* 1 means to print initial x and d
859 c$$$c iv(24)=1 ! DEBUG
860 c$$$* min val for v(radfac) default is 0.1
862 c$$$* max val for v(radfac) default is 4.0
865 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
866 c$$$* the sumsl default is 0.1
868 c$$$* false conv if (act fnctn decrease) .lt. v(34)
869 c$$$* the sumsl default is 100*machep
870 c$$$ v(34)=v(34)/100.0D0
871 c$$$* absolute convergence
872 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
874 c$$$* relative convergence
875 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
877 c$$$* controls initial step size
879 c$$$* large vals of d correspond to small components of step
883 c$$$ do i=nphi+1,nvar
887 c$$$ call geom_to_var(nvar,x)
888 c$$$ IF (mask_r) THEN
889 c$$$ do i=1,nres ! Just in case...
893 c$$$ call x2xx(x,xx,nvar_restr)
894 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
895 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
898 c$$$c When minimizing ALL side-chains, etotal_sc is a little
899 c$$$c faster if we don't set mask_r
905 c$$$ call x2xx(x,xx,nvar_restr)
906 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
907 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
910 c$$$ call var_to_geom(nvar,x)
911 c$$$ call chainbuild_sc
918 c$$$C--------------------------------------------------------------------------
920 c$$$ subroutine chainbuild_sc
922 c$$$ include 'DIMENSIONS'
923 c$$$ include 'COMMON.VAR'
924 c$$$ include 'COMMON.INTERACT'
926 c$$$c Local variables
931 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
932 c$$$ call locate_side_chain(i)
939 c$$$C--------------------------------------------------------------------------
941 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
945 c$$$ include 'DIMENSIONS'
946 c$$$ include 'COMMON.DERIV'
947 c$$$ include 'COMMON.VAR'
948 c$$$ include 'COMMON.MINIM'
949 c$$$ include 'COMMON.IOUNITS'
951 c$$$c Input arguments
953 c$$$ double precision x(maxvar)
954 c$$$ double precision ufparm
957 c$$$c Input/Output arguments
959 c$$$ integer uiparm(1)
960 c$$$ double precision urparm(1)
962 c$$$c Output arguments
963 c$$$ double precision f
965 c$$$c Local variables
966 c$$$ double precision energia(0:n_ene)
968 c$$$c Variables used to intercept NaNs
969 c$$$ double precision x_sum
978 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
981 c$$$ x_sum=x_sum+x(i_NAN)
983 c$$$c Calculate the energy only if the coordinates are ok
984 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
985 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
991 c$$$ call var_to_geom_restr(n,x)
993 c$$$ call chainbuild_sc
994 c$$$ call etotal_sc(energia(0))
996 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
1005 c$$$c-------------------------------------------------------
1007 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
1011 c$$$ include 'DIMENSIONS'
1012 c$$$ include 'COMMON.CHAIN'
1013 c$$$ include 'COMMON.DERIV'
1014 c$$$ include 'COMMON.VAR'
1015 c$$$ include 'COMMON.INTERACT'
1016 c$$$ include 'COMMON.MINIM'
1018 c$$$c Input arguments
1020 c$$$ double precision x(maxvar)
1021 c$$$ double precision ufparm
1022 c$$$ external ufparm
1024 c$$$c Input/Output arguments
1026 c$$$ integer uiparm(1)
1027 c$$$ double precision urparm(1)
1029 c$$$c Output arguments
1030 c$$$ double precision g(maxvar)
1032 c$$$c Local variables
1033 c$$$ double precision f,gphii,gthetai,galphai,gomegai
1034 c$$$ integer ig,ind,i,j,k,igall,ij
1037 c$$$ icg=mod(nf,2)+1
1038 c$$$ if (nf-nfl+1) 20,30,40
1039 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1040 c$$$c write (iout,*) 'grad 20'
1041 c$$$ if (nf.eq.0) return
1043 c$$$ 30 call var_to_geom_restr(n,x)
1044 c$$$ call chainbuild_sc
1046 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1048 c$$$ 40 call cartder
1050 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1056 c$$$ IF (mask_phi(i+2).eq.1) THEN
1058 c$$$ do j=i+1,nres-1
1061 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1062 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1068 c$$$ ind=ind+nres-1-i
1075 c$$$ IF (mask_theta(i+2).eq.1) THEN
1078 c$$$ do j=i+1,nres-1
1081 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1082 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1087 c$$$ ind=ind+nres-1-i
1092 c$$$ if (itype(i).ne.10) then
1093 c$$$ IF (mask_side(i).eq.1) THEN
1097 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1106 c$$$ if (itype(i).ne.10) then
1107 c$$$ IF (mask_side(i).eq.1) THEN
1111 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1119 c$$$C Add the components corresponding to local energy terms.
1126 c$$$ if (mask_phi(i).eq.1) then
1128 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1134 c$$$ if (mask_theta(i).eq.1) then
1136 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1142 c$$$ if (itype(i).ne.10) then
1144 c$$$ if (mask_side(i).eq.1) then
1146 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1153 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1159 c$$$C-----------------------------------------------------------------------------
1161 c$$$ subroutine etotal_sc(energy_sc)
1165 c$$$ include 'DIMENSIONS'
1166 c$$$ include 'COMMON.VAR'
1167 c$$$ include 'COMMON.INTERACT'
1168 c$$$ include 'COMMON.DERIV'
1169 c$$$ include 'COMMON.FFIELD'
1171 c$$$c Output arguments
1172 c$$$ double precision energy_sc(0:n_ene)
1174 c$$$c Local variables
1175 c$$$ double precision evdw,escloc
1180 c$$$ energy_sc(i)=0.0D0
1183 c$$$ if (mask_r) then
1184 c$$$ call egb_sc(evdw)
1185 c$$$ call esc_sc(escloc)
1188 c$$$ call esc(escloc)
1191 c$$$ if (evdw.eq.1.0D20) then
1192 c$$$ energy_sc(0)=evdw
1194 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1196 c$$$ energy_sc(1)=evdw
1197 c$$$ energy_sc(12)=escloc
1200 c$$$C Sum up the components of the Cartesian gradient.
1204 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1211 c$$$C-----------------------------------------------------------------------------
1213 c$$$ subroutine egb_sc(evdw)
1215 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1216 c$$$C assuming the Gay-Berne potential of interaction.
1218 c$$$ implicit real*8 (a-h,o-z)
1219 c$$$ include 'DIMENSIONS'
1220 c$$$ include 'COMMON.GEO'
1221 c$$$ include 'COMMON.VAR'
1222 c$$$ include 'COMMON.LOCAL'
1223 c$$$ include 'COMMON.CHAIN'
1224 c$$$ include 'COMMON.DERIV'
1225 c$$$ include 'COMMON.NAMES'
1226 c$$$ include 'COMMON.INTERACT'
1227 c$$$ include 'COMMON.IOUNITS'
1228 c$$$ include 'COMMON.CALC'
1229 c$$$ include 'COMMON.CONTROL'
1232 c$$$ energy_dec=.false.
1233 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1236 c$$$c if (icall.eq.0) lprn=.false.
1238 c$$$ do i=iatsc_s,iatsc_e
1240 c$$$ itypi1=itype(i+1)
1244 c$$$ dxi=dc_norm(1,nres+i)
1245 c$$$ dyi=dc_norm(2,nres+i)
1246 c$$$ dzi=dc_norm(3,nres+i)
1247 c$$$c dsci_inv=dsc_inv(itypi)
1248 c$$$ dsci_inv=vbld_inv(i+nres)
1249 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1250 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1252 c$$$C Calculate SC interaction energy.
1254 c$$$ do iint=1,nint_gr(i)
1255 c$$$ do j=istart(i,iint),iend(i,iint)
1256 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1259 c$$$c dscj_inv=dsc_inv(itypj)
1260 c$$$ dscj_inv=vbld_inv(j+nres)
1261 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1262 c$$$c & 1.0d0/vbld(j+nres)
1263 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1264 c$$$ sig0ij=sigma(itypi,itypj)
1265 c$$$ chi1=chi(itypi,itypj)
1266 c$$$ chi2=chi(itypj,itypi)
1267 c$$$ chi12=chi1*chi2
1268 c$$$ chip1=chip(itypi)
1269 c$$$ chip2=chip(itypj)
1270 c$$$ chip12=chip1*chip2
1271 c$$$ alf1=alp(itypi)
1272 c$$$ alf2=alp(itypj)
1273 c$$$ alf12=0.5D0*(alf1+alf2)
1274 c$$$C For diagnostics only!!!
1284 c$$$ xj=c(1,nres+j)-xi
1285 c$$$ yj=c(2,nres+j)-yi
1286 c$$$ zj=c(3,nres+j)-zi
1287 c$$$ dxj=dc_norm(1,nres+j)
1288 c$$$ dyj=dc_norm(2,nres+j)
1289 c$$$ dzj=dc_norm(3,nres+j)
1290 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1291 c$$$c write (iout,*) "j",j," dc_norm",
1292 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1293 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1294 c$$$ rij=dsqrt(rrij)
1295 c$$$C Calculate angle-dependent terms of energy and contributions to their
1297 c$$$ call sc_angular
1298 c$$$ sigsq=1.0D0/sigsq
1299 c$$$ sig=sig0ij*dsqrt(sigsq)
1300 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1301 c$$$c for diagnostics; uncomment
1302 c$$$c rij_shift=1.2*sig0ij
1303 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1304 c$$$ if (rij_shift.le.0.0D0) then
1306 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1307 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1308 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1311 c$$$ sigder=-sig*sigsq
1312 c$$$c---------------------------------------------------------------
1313 c$$$ rij_shift=1.0D0/rij_shift
1314 c$$$ fac=rij_shift**expon
1315 c$$$ e1=fac*fac*aa(itypi,itypj)
1316 c$$$ e2=fac*bb(itypi,itypj)
1317 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1318 c$$$ eps2der=evdwij*eps3rt
1319 c$$$ eps3der=evdwij*eps2rt
1320 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1321 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1322 c$$$ evdwij=evdwij*eps2rt*eps3rt
1323 c$$$ evdw=evdw+evdwij
1325 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1326 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1327 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1328 c$$$ & restyp(itypi),i,restyp(itypj),j,
1329 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1330 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1331 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1335 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1336 c$$$ & 'evdw',i,j,evdwij
1338 c$$$C Calculate gradient components.
1339 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1340 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1341 c$$$ sigder=fac*sigder
1344 c$$$C Calculate the radial part of the gradient
1348 c$$$C Calculate angular part of the gradient.
1354 c$$$ energy_dec=.false.
1358 c$$$c-----------------------------------------------------------------------------
1360 c$$$ subroutine esc_sc(escloc)
1361 c$$$C Calculate the local energy of a side chain and its derivatives in the
1362 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1363 c$$$C ALPHA and OMEGA.
1364 c$$$ implicit real*8 (a-h,o-z)
1365 c$$$ include 'DIMENSIONS'
1366 c$$$ include 'COMMON.GEO'
1367 c$$$ include 'COMMON.LOCAL'
1368 c$$$ include 'COMMON.VAR'
1369 c$$$ include 'COMMON.INTERACT'
1370 c$$$ include 'COMMON.DERIV'
1371 c$$$ include 'COMMON.CHAIN'
1372 c$$$ include 'COMMON.IOUNITS'
1373 c$$$ include 'COMMON.NAMES'
1374 c$$$ include 'COMMON.FFIELD'
1375 c$$$ include 'COMMON.CONTROL'
1376 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1377 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1378 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1379 c$$$ delta=0.02d0*pi
1381 c$$$c write (iout,'(a)') 'ESC'
1382 c$$$ do i=loc_start,loc_end
1383 c$$$ IF (mask_side(i).eq.1) THEN
1385 c$$$ if (it.eq.10) goto 1
1386 c$$$ nlobit=nlob(it)
1387 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1388 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1389 c$$$ theti=theta(i+1)-pipol
1390 c$$$ x(1)=dtan(theti)
1394 c$$$ if (x(2).gt.pi-delta) then
1396 c$$$ xtemp(2)=pi-delta
1398 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1400 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1401 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1402 c$$$ & escloci,dersc(2))
1403 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1404 c$$$ & ddersc0(1),dersc(1))
1405 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1406 c$$$ & ddersc0(3),dersc(3))
1407 c$$$ xtemp(2)=pi-delta
1408 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1410 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1411 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1412 c$$$ & dersc0(2),esclocbi,dersc02)
1413 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1414 c$$$ & dersc12,dersc01)
1415 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1416 c$$$ dersc0(1)=dersc01
1417 c$$$ dersc0(2)=dersc02
1418 c$$$ dersc0(3)=0.0d0
1420 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1422 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1423 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1424 c$$$c & esclocbi,ss,ssd
1425 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1426 c$$$c escloci=esclocbi
1427 c$$$c write (iout,*) escloci
1428 c$$$ else if (x(2).lt.delta) then
1432 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1434 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1435 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1436 c$$$ & escloci,dersc(2))
1437 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1438 c$$$ & ddersc0(1),dersc(1))
1439 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1440 c$$$ & ddersc0(3),dersc(3))
1442 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1444 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1445 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1446 c$$$ & dersc0(2),esclocbi,dersc02)
1447 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1448 c$$$ & dersc12,dersc01)
1449 c$$$ dersc0(1)=dersc01
1450 c$$$ dersc0(2)=dersc02
1451 c$$$ dersc0(3)=0.0d0
1452 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1454 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1456 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1457 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1458 c$$$c & esclocbi,ss,ssd
1459 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1460 c$$$c write (iout,*) escloci
1462 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1465 c$$$ escloc=escloc+escloci
1466 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1467 c$$$ & 'escloc',i,escloci
1468 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1470 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1471 c$$$ & wscloc*dersc(1)
1472 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1473 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1480 c$$$C-----------------------------------------------------------------------------
1482 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1484 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1485 c$$$C assuming the Gay-Berne potential of interaction.
1487 c$$$ implicit real*8 (a-h,o-z)
1488 c$$$ include 'DIMENSIONS'
1489 c$$$ include 'COMMON.GEO'
1490 c$$$ include 'COMMON.VAR'
1491 c$$$ include 'COMMON.LOCAL'
1492 c$$$ include 'COMMON.CHAIN'
1493 c$$$ include 'COMMON.DERIV'
1494 c$$$ include 'COMMON.NAMES'
1495 c$$$ include 'COMMON.INTERACT'
1496 c$$$ include 'COMMON.IOUNITS'
1497 c$$$ include 'COMMON.CALC'
1498 c$$$ include 'COMMON.CONTROL'
1501 c$$$ energy_dec=.false.
1502 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1506 c$$$c$$$ do i=iatsc_s,iatsc_e
1509 c$$$ itypi1=itype(i+1)
1513 c$$$ dxi=dc_norm(1,nres+i)
1514 c$$$ dyi=dc_norm(2,nres+i)
1515 c$$$ dzi=dc_norm(3,nres+i)
1516 c$$$c dsci_inv=dsc_inv(itypi)
1517 c$$$ dsci_inv=vbld_inv(i+nres)
1518 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1519 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1521 c$$$C Calculate SC interaction energy.
1523 c$$$c$$$ do iint=1,nint_gr(i)
1524 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1528 c$$$c dscj_inv=dsc_inv(itypj)
1529 c$$$ dscj_inv=vbld_inv(j+nres)
1530 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1531 c$$$c & 1.0d0/vbld(j+nres)
1532 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1533 c$$$ sig0ij=sigma(itypi,itypj)
1534 c$$$ chi1=chi(itypi,itypj)
1535 c$$$ chi2=chi(itypj,itypi)
1536 c$$$ chi12=chi1*chi2
1537 c$$$ chip1=chip(itypi)
1538 c$$$ chip2=chip(itypj)
1539 c$$$ chip12=chip1*chip2
1540 c$$$ alf1=alp(itypi)
1541 c$$$ alf2=alp(itypj)
1542 c$$$ alf12=0.5D0*(alf1+alf2)
1543 c$$$C For diagnostics only!!!
1553 c$$$ xj=c(1,nres+j)-xi
1554 c$$$ yj=c(2,nres+j)-yi
1555 c$$$ zj=c(3,nres+j)-zi
1556 c$$$ dxj=dc_norm(1,nres+j)
1557 c$$$ dyj=dc_norm(2,nres+j)
1558 c$$$ dzj=dc_norm(3,nres+j)
1559 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1560 c$$$c write (iout,*) "j",j," dc_norm",
1561 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1562 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1563 c$$$ rij=dsqrt(rrij)
1564 c$$$C Calculate angle-dependent terms of energy and contributions to their
1566 c$$$ call sc_angular
1567 c$$$ sigsq=1.0D0/sigsq
1568 c$$$ sig=sig0ij*dsqrt(sigsq)
1569 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1570 c$$$c for diagnostics; uncomment
1571 c$$$c rij_shift=1.2*sig0ij
1572 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1573 c$$$ if (rij_shift.le.0.0D0) then
1575 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1576 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1577 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1580 c$$$ sigder=-sig*sigsq
1581 c$$$c---------------------------------------------------------------
1582 c$$$ rij_shift=1.0D0/rij_shift
1583 c$$$ fac=rij_shift**expon
1584 c$$$ e1=fac*fac*aa(itypi,itypj)
1585 c$$$ e2=fac*bb(itypi,itypj)
1586 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1587 c$$$ eps2der=evdwij*eps3rt
1588 c$$$ eps3der=evdwij*eps2rt
1589 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1590 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1591 c$$$ evdwij=evdwij*eps2rt*eps3rt
1592 c$$$ evdw=evdw+evdwij
1594 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1595 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1596 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1597 c$$$ & restyp(itypi),i,restyp(itypj),j,
1598 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1599 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1600 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1604 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1605 c$$$ & 'evdw',i,j,evdwij
1607 c$$$C Calculate gradient components.
1608 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1609 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1610 c$$$ sigder=fac*sigder
1613 c$$$C Calculate the radial part of the gradient
1617 c$$$C Calculate angular part of the gradient.
1620 c$$$c$$$ enddo ! iint
1622 c$$$ energy_dec=.false.
1626 c$$$C-----------------------------------------------------------------------------
1628 c$$$ subroutine perturb_side_chain(i,angle)
1632 c$$$ include 'DIMENSIONS'
1633 c$$$ include 'COMMON.CHAIN'
1634 c$$$ include 'COMMON.GEO'
1635 c$$$ include 'COMMON.VAR'
1636 c$$$ include 'COMMON.LOCAL'
1637 c$$$ include 'COMMON.IOUNITS'
1639 c$$$c External functions
1640 c$$$ external ran_number
1641 c$$$ double precision ran_number
1643 c$$$c Input arguments
1645 c$$$ double precision angle ! In degrees
1647 c$$$c Local variables
1649 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1653 c$$$ rad_ang=angle*deg2rad
1656 c$$$ do while (length.lt.0.01)
1657 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1658 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1659 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1660 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1661 c$$$ + rand_v(3)*rand_v(3)
1662 c$$$ length=sqrt(length)
1663 c$$$ rand_v(1)=rand_v(1)/length
1664 c$$$ rand_v(2)=rand_v(2)/length
1665 c$$$ rand_v(3)=rand_v(3)/length
1666 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1667 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1668 c$$$ length=1.0D0-cost*cost
1669 c$$$ if (length.lt.0.0D0) length=0.0D0
1670 c$$$ length=sqrt(length)
1671 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1672 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1673 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1675 c$$$ rand_v(1)=rand_v(1)/length
1676 c$$$ rand_v(2)=rand_v(2)/length
1677 c$$$ rand_v(3)=rand_v(3)/length
1679 c$$$ cost=dcos(rad_ang)
1680 c$$$ sint=dsin(rad_ang)
1681 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1682 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1683 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1684 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1685 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1686 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1687 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1688 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1689 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1691 c$$$ call chainbuild_cart
1696 c$$$c----------------------------------------------------------------------------
1698 c$$$ subroutine ss_relax3(i_in,j_in)
1702 c$$$ include 'DIMENSIONS'
1703 c$$$ include 'COMMON.VAR'
1704 c$$$ include 'COMMON.CHAIN'
1705 c$$$ include 'COMMON.IOUNITS'
1706 c$$$ include 'COMMON.INTERACT'
1708 c$$$c External functions
1709 c$$$ external ran_number
1710 c$$$ double precision ran_number
1712 c$$$c Input arguments
1713 c$$$ integer i_in,j_in
1715 c$$$c Local variables
1716 c$$$ double precision energy_sc(0:n_ene),etot
1717 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1718 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1719 c$$$ integer n,i_pert,i
1720 c$$$ logical notdone
1729 c$$$ mask_side(i_in)=1
1730 c$$$ mask_side(j_in)=1
1732 c$$$ call etotal_sc(energy_sc)
1733 c$$$ etot=energy_sc(0)
1734 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1735 c$$$c + energy_sc(1),energy_sc(12)
1739 c$$$ do while (notdone)
1740 c$$$ if (mod(n,2).eq.0) then
1748 c$$$ org_dc(i)=dc(i,i_pert+nres)
1749 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1750 c$$$ org_c(i)=c(i,i_pert+nres)
1752 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1753 c$$$ call perturb_side_chain(i_pert,ang_pert)
1754 c$$$ call etotal_sc(energy_sc)
1755 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1756 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1757 c$$$ if (rand_fact.lt.exp_fact) then
1758 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1759 c$$$c + energy_sc(1),energy_sc(12)
1760 c$$$ etot=energy_sc(0)
1762 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1763 c$$$c + energy_sc(1),energy_sc(12)
1765 c$$$ dc(i,i_pert+nres)=org_dc(i)
1766 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1767 c$$$ c(i,i_pert+nres)=org_c(i)
1771 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1779 c$$$c----------------------------------------------------------------------------
1781 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1783 c$$$ include 'DIMENSIONS'
1785 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1786 c$$$*********************************************************************
1787 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1788 c$$$* the calling subprogram. *
1789 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1790 c$$$* calculated in the usual pythagorean way. *
1791 c$$$* absolute convergence occurs when the function is within v(31) of *
1792 c$$$* zero. unless you know the minimum value in advance, abs convg *
1793 c$$$* is probably not useful. *
1794 c$$$* relative convergence is when the model predicts that the function *
1795 c$$$* will decrease by less than v(32)*abs(fun). *
1796 c$$$*********************************************************************
1797 c$$$ include 'COMMON.IOUNITS'
1798 c$$$ include 'COMMON.VAR'
1799 c$$$ include 'COMMON.GEO'
1800 c$$$ include 'COMMON.MINIM'
1801 c$$$ include 'COMMON.CHAIN'
1803 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1804 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1805 c$$$ + orig_ss_dist(maxres2,maxres2)
1807 c$$$ double precision etot
1808 c$$$ integer iretcode,nfun,i_in,j_in
1811 c$$$ double precision dist
1812 c$$$ external ss_func,fdum
1813 c$$$ double precision ss_func,fdum
1815 c$$$ integer iv(liv),uiparm(2)
1816 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1820 c$$$ call deflt(2,iv,liv,lv,v)
1821 c$$$* 12 means fresh start, dont call deflt
1823 c$$$* max num of fun calls
1824 c$$$ if (maxfun.eq.0) maxfun=500
1826 c$$$* max num of iterations
1827 c$$$ if (maxmin.eq.0) maxmin=1000
1829 c$$$* controls output
1831 c$$$* selects output unit
1834 c$$$* 1 means to print out result
1836 c$$$* 1 means to print out summary stats
1838 c$$$* 1 means to print initial x and d
1840 c$$$* min val for v(radfac) default is 0.1
1842 c$$$* max val for v(radfac) default is 4.0
1845 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1846 c$$$* the sumsl default is 0.1
1848 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1849 c$$$* the sumsl default is 100*machep
1850 c$$$ v(34)=v(34)/100.0D0
1851 c$$$* absolute convergence
1852 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1855 c$$$* relative convergence
1856 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1859 c$$$* controls initial step size
1861 c$$$* large vals of d correspond to small components of step
1868 c$$$ orig_ss_dc(j,i)=dc(j,i)
1871 c$$$ call geom_to_var(nvar,orig_ss_var)
1875 c$$$ orig_ss_dist(j,i)=dist(j,i)
1876 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1877 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1878 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1890 c$$$ if (ialph(i,1).gt.0) then
1893 c$$$ x(k)=dc(j,i+nres)
1900 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1903 c$$$ nfun=iv(6)+iv(30)
1913 c$$$ if (ialph(i,1).gt.0) then
1916 c$$$ dc(j,i+nres)=x(k)
1920 c$$$ call chainbuild_cart
1925 c$$$C-----------------------------------------------------------------------------
1927 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1929 c$$$ include 'DIMENSIONS'
1930 c$$$ include 'COMMON.DERIV'
1931 c$$$ include 'COMMON.IOUNITS'
1932 c$$$ include 'COMMON.VAR'
1933 c$$$ include 'COMMON.CHAIN'
1934 c$$$ include 'COMMON.INTERACT'
1935 c$$$ include 'COMMON.SBRIDGE'
1937 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1938 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1939 c$$$ + orig_ss_dist(maxres2,maxres2)
1942 c$$$ double precision x(maxres6)
1944 c$$$ double precision f
1945 c$$$ integer uiparm(2)
1946 c$$$ real*8 urparm(1)
1947 c$$$ external ufparm
1948 c$$$ double precision ufparm
1951 c$$$ double precision dist
1953 c$$$ integer i,j,k,ss_i,ss_j
1954 c$$$ double precision tempf,var(maxvar)
1969 c$$$ if (ialph(i,1).gt.0) then
1972 c$$$ dc(j,i+nres)=x(k)
1976 c$$$ call chainbuild_cart
1978 c$$$ call geom_to_var(nvar,var)
1980 c$$$c Constraints on all angles
1982 c$$$ tempf=var(i)-orig_ss_var(i)
1983 c$$$ f=f+tempf*tempf
1986 c$$$c Constraints on all distances
1988 c$$$ if (i.gt.1) then
1989 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1990 c$$$ f=f+tempf*tempf
1993 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1994 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1995 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1996 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1997 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1998 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1999 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
2000 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2004 c$$$c Constraints for the relevant CYS-CYS
2005 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
2006 c$$$ f=f+tempf*tempf
2007 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
2009 c$$$c$$$ if (nf.ne.nfl) then
2010 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
2011 c$$$c$$$ + f,dist(5+nres,14+nres)
2019 c$$$C-----------------------------------------------------------------------------
2020 c$$$C-----------------------------------------------------------------------------
2021 subroutine triple_ssbond_ene(resi,resj,resk,eij)
2022 include 'DIMENSIONS'
2023 include 'COMMON.SBRIDGE'
2024 include 'COMMON.CHAIN'
2025 include 'COMMON.DERIV'
2026 include 'COMMON.LOCAL'
2027 include 'COMMON.INTERACT'
2028 include 'COMMON.VAR'
2029 include 'COMMON.IOUNITS'
2030 include 'COMMON.CALC'
2033 C include 'COMMON.MD'
2037 c External functions
2038 double precision h_base
2042 integer resi,resj,resk
2045 double precision eij,eij1,eij2,eij3
2049 c integer itypi,itypj,k,l
2050 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2051 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2052 double precision xik,yik,zik,xjk,yjk,zjk
2053 double precision sig0ij,ljd,sig,fac,e1,e2
2054 double precision dcosom1(3),dcosom2(3),ed
2055 double precision pom1,pom2
2056 double precision ljA,ljB,ljXs
2057 double precision d_ljB(1:3)
2058 double precision ssA,ssB,ssC,ssXs
2059 double precision ssxm,ljxm,ssm,ljm
2060 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2065 C write(iout,*) resi,resj,resk
2067 dxi=dc_norm(1,nres+i)
2068 dyi=dc_norm(2,nres+i)
2069 dzi=dc_norm(3,nres+i)
2070 dsci_inv=vbld_inv(i+nres)
2080 dxj=dc_norm(1,nres+j)
2081 dyj=dc_norm(2,nres+j)
2082 dzj=dc_norm(3,nres+j)
2083 dscj_inv=vbld_inv(j+nres)
2089 dxk=dc_norm(1,nres+k)
2090 dyk=dc_norm(2,nres+k)
2091 dzk=dc_norm(3,nres+k)
2092 dscj_inv=vbld_inv(k+nres)
2102 rrij=(xij*xij+yij*yij+zij*zij)
2103 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2104 rrik=(xik*xik+yik*yik+zik*zik)
2106 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2108 C there are three combination of distances for each trisulfide bonds
2109 C The first case the ith atom is the center
2110 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2111 C distance y is second distance the a,b,c,d are parameters derived for
2112 C this problem d parameter was set as a penalty currenlty set to 1.
2113 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2114 C second case jth atom is center
2115 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2116 C the third case kth atom is the center
2117 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2122 C write(iout,*)i,j,k,eij
2123 C The energy penalty calculated now time for the gradient part
2124 C derivative over rij
2125 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2126 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2131 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2132 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2135 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2136 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2138 C now derivative over rik
2139 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2140 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2145 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2146 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2149 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2150 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2152 C now derivative over rjk
2153 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2154 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2159 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2160 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2163 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2164 gvdwc(l,k)=gvdwc(l,k)+gg(l)