1 c----------------------------------------------------------------------------
2 subroutine check_energies
6 implicit real*8 (a-h,o-z)
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)
91 implicit real*8 (a-h,o-z)
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 integer xshift,yshift,zshift
134 c-------END FIRST METHOD
135 c-------SECOND METHOD
136 c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
137 c-------END SECOND METHOD
140 logical checkstop,transgrad
141 common /sschecks/ checkstop,transgrad
143 integer icheck,nicheck,jcheck,njcheck
144 double precision echeck(-1:1),deps,ssx0,ljx0
145 c-------END TESTING CODE
152 dxi=dc_norm(1,nres+i)
153 dyi=dc_norm(2,nres+i)
154 dzi=dc_norm(3,nres+i)
155 dsci_inv=vbld_inv(i+nres)
160 if (xi.lt.0) xi=xi+boxxsize
162 if (yi.lt.0) yi=yi+boxysize
164 if (zi.lt.0) zi=zi+boxzsize
165 C define scaling factor for lipids
167 C if (positi.le.0) positi=positi+boxzsize
169 C first for peptide groups
170 c for each residue check if it is in lipid or lipid water border area
171 if ((zi.gt.bordlipbot)
172 &.and.(zi.lt.bordliptop)) then
173 C the energy transfer exist
174 if (zi.lt.buflipbot) then
175 C what fraction I am in
177 & ((positi-bordlipbot)/lipbufthick)
178 C lipbufthick is thickenes of lipid buffore
179 sslipi=sscalelip(fracinbuf)
180 ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
181 elseif (zi.gt.bufliptop) then
182 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
183 sslipi=sscalelip(fracinbuf)
184 ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
199 if (xj.lt.0) xj=xj+boxxsize
201 if (yj.lt.0) yj=yj+boxysize
203 if (zj.lt.0) zj=zj+boxzsize
204 if ((zj.gt.bordlipbot)
205 &.and.(zj.lt.bordliptop)) then
206 C the energy transfer exist
207 if (zj.lt.buflipbot) then
208 C what fraction I am in
210 & ((positi-bordlipbot)/lipbufthick)
211 C lipbufthick is thickenes of lipid buffore
212 sslipj=sscalelip(fracinbuf)
213 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
214 elseif (zi.gt.bufliptop) then
215 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
216 sslipj=sscalelip(fracinbuf)
217 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
226 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
227 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
228 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
229 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
231 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
243 xj=xj_safe+xshift*boxxsize
244 yj=yj_safe+yshift*boxysize
245 zj=zj_safe+zshift*boxzsize
246 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
247 if(dist_temp.lt.dist_init) then
257 if (subchap.eq.1) then
267 dxj=dc_norm(1,nres+j)
268 dyj=dc_norm(2,nres+j)
269 dzj=dc_norm(3,nres+j)
270 dscj_inv=vbld_inv(j+nres)
272 chi1=chi(itypi,itypj)
273 chi2=chi(itypj,itypi)
280 alf12=0.5D0*(alf1+alf2)
282 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
283 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
284 c The following are set in sc_angular
288 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
289 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
290 c om12=dxi*dxj+dyi*dyj+dzi*dzj
292 rij=1.0D0/rij ! Reset this so it makes sense
294 sig0ij=sigma(itypi,itypj)
295 sig=sig0ij*dsqrt(1.0D0/sigsq)
298 ljA=eps1*eps2rt**2*eps3rt**2
301 ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
306 deltat12=om2-om1+2.0d0
311 & +akth*(deltat1*deltat1+deltat2*deltat2)
312 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
313 ssxm=ssXs-0.5D0*ssB/ssA
316 c$$$c Some extra output
317 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
318 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
319 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
320 c$$$ if (ssx0.gt.0.0d0) then
321 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
325 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
326 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
327 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
329 c-------END TESTING CODE
332 c Stop and plot energy and derivative as a function of distance
334 ssm=ssC-0.25D0*ssB*ssB/ssA
335 ljm=-0.25D0*ljB*bb/aa
337 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
345 if (.not.checkstop) then
352 if (checkstop) rij=(ssxm-1.0d0)+
353 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
354 c-------END TESTING CODE
356 if (rij.gt.ljxm) then
359 fac=(1.0D0/ljd)**expon
362 eij=eps1*eps2rt*eps3rt*(e1+e2)
363 C write(iout,*) eij,'TU?1'
366 eij=eij*eps2rt*eps3rt
369 e1=e1*eps1*eps2rt**2*eps3rt**2
370 ed=-expon*(e1+eij)/ljd
372 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
373 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
374 eom12=eij*eps1_om12+eps2der*eps2rt_om12
375 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
376 else if (rij.lt.ssxm) then
379 eij=ssA*ssd*ssd+ssB*ssd+ssC
380 C write(iout,*) 'TU?2',ssc,ssd
381 ed=2*akcm*ssd+akct*deltat12
383 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
384 eom1=-2*akth*deltat1-pom1-om2*pom2
385 eom2= 2*akth*deltat2+pom1-om1*pom2
388 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
390 d_ssxm(1)=0.5D0*akct/ssA
394 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
395 d_ljxm(2)=d_ljxm(1)*sigsq_om2
396 d_ljxm(3)=d_ljxm(1)*sigsq_om12
397 d_ljxm(1)=d_ljxm(1)*sigsq_om1
399 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
402 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
406 ssm=ssC-0.25D0*ssB*ssB/ssA
407 d_ssm(1)=0.5D0*akct*ssB/ssA
408 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
409 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
411 f1=(rij-xm)/(ssxm-xm)
412 f2=(rij-ssxm)/(xm-ssxm)
416 C write(iout,*) eij,'TU?3'
417 delta_inv=1.0d0/(xm-ssxm)
418 deltasq_inv=delta_inv*delta_inv
420 fac1=deltasq_inv*fac*(xm-rij)
421 fac2=deltasq_inv*fac*(rij-ssxm)
422 ed=delta_inv*(Ht*hd2-ssm*hd1)
423 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
424 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
425 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
428 ljm=-0.25D0*ljB*bb/aa
429 d_ljm(1)=-0.5D0*bb/aa*ljB
430 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
431 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
433 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
434 f1=(rij-ljxm)/(xm-ljxm)
435 f2=(rij-xm)/(ljxm-xm)
439 C write(iout,*) 'TU?4',ssA
440 delta_inv=1.0d0/(ljxm-xm)
441 deltasq_inv=delta_inv*delta_inv
443 fac1=deltasq_inv*fac*(ljxm-rij)
444 fac2=deltasq_inv*fac*(rij-xm)
445 ed=delta_inv*(ljm*hd2-Ht*hd1)
446 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
447 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
448 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
450 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
452 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
458 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
459 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
460 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
462 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
463 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
464 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
465 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
468 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
470 c$$$ d_ljm(k)=ljm*d_ljB(k)
474 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
475 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
476 c$$$ d_ss(2)=akct*ssd
477 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
478 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
481 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
482 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
483 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
485 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
486 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
488 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
490 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
491 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
492 c$$$ h1=h_base(f1,hd1)
493 c$$$ h2=h_base(f2,hd2)
494 c$$$ eij=ss*h1+ljf*h2
495 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
496 c$$$ deltasq_inv=delta_inv*delta_inv
497 c$$$ fac=ljf*hd2-ss*hd1
498 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
499 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
500 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
501 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
502 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
503 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
504 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
506 c$$$ havebond=.false.
507 c$$$ if (ed.gt.0.0d0) havebond=.true.
508 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
511 C write(iout,*) 'havebond',havebond
515 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
516 c write(iout,'(a15,f12.2,f8.1,2i5)')
517 c & "SSBOND_E_FORM",totT,t_bath,i,j
521 dyn_ssbond_ij(i,j)=eij
522 else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
523 dyn_ssbond_ij(i,j)=1.0d300
526 c write(iout,'(a15,f12.2,f8.1,2i5)')
527 c & "SSBOND_E_BREAK",totT,t_bath,i,j
534 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
535 & "CHECKSTOP",rij,eij,ed
540 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
547 c-------END TESTING CODE
550 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
551 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
554 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
557 gvdwx(k,i)=gvdwx(k,i)-gg(k)
558 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
559 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
560 gvdwx(k,j)=gvdwx(k,j)+gg(k)
561 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
562 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
566 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
571 gvdwc(l,i)=gvdwc(l,i)-gg(l)
572 gvdwc(l,j)=gvdwc(l,j)+gg(l)
578 C-----------------------------------------------------------------------------
580 double precision function h_base(x,deriv)
581 c A smooth function going 0->1 in range [0,1]
582 c It should NOT be called outside range [0,1], it will not work there.
589 double precision deriv
595 c Two parabolas put together. First derivative zero at extrema
596 c$$$ if (x.lt.0.5D0) then
597 c$$$ h_base=2.0D0*x*x
601 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
602 c$$$ deriv=4.0D0*deriv
605 c Third degree polynomial. First derivative zero at extrema
606 h_base=x*x*(3.0d0-2.0d0*x)
607 deriv=6.0d0*x*(1.0d0-x)
609 c Fifth degree polynomial. First and second derivatives zero at extrema
611 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
613 c$$$ deriv=deriv*deriv
614 c$$$ deriv=30.0d0*xsq*deriv
619 c----------------------------------------------------------------------------
621 subroutine dyn_set_nss
622 c Adjust nss and other relevant variables based on dyn_ssbond_ij
626 implicit real*8 (a-h,o-z)
631 include 'COMMON.SBRIDGE'
632 include 'COMMON.CHAIN'
633 include 'COMMON.IOUNITS'
634 C include 'COMMON.SETUP'
637 C include 'COMMON.MD'
642 double precision emin
644 integer diff,allflag(maxdim),allnss,
645 & allihpb(maxdim),alljhpb(maxdim),
646 & newnss,newihpb(maxdim),newjhpb(maxdim)
648 integer i_newnss(1024),displ(0:1024)
649 integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
654 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
663 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
667 if (allflag(i).eq.0 .and.
668 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
669 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
673 if (emin.lt.1.0d300) then
676 if (allflag(i).eq.0 .and.
677 & (allihpb(i).eq.allihpb(imin) .or.
678 & alljhpb(i).eq.allihpb(imin) .or.
679 & allihpb(i).eq.alljhpb(imin) .or.
680 & alljhpb(i).eq.alljhpb(imin))) then
687 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
691 if (allflag(i).eq.1) then
693 newihpb(newnss)=allihpb(i)
694 newjhpb(newnss)=alljhpb(i)
699 if (nfgtasks.gt.1)then
701 call MPI_Reduce(newnss,g_newnss,1,
702 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
703 call MPI_Gather(newnss,1,MPI_INTEGER,
704 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
707 displ(i)=i_newnss(i-1)+displ(i-1)
709 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
710 & g_newihpb,i_newnss,displ,MPI_INTEGER,
712 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
713 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
715 if(fg_rank.eq.0) then
716 c print *,'g_newnss',g_newnss
717 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
718 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
721 newihpb(i)=g_newihpb(i)
722 newjhpb(i)=g_newjhpb(i)
730 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
735 if (idssb(i).eq.newihpb(j) .and.
736 & jdssb(i).eq.newjhpb(j)) found=.true.
740 c if (.not.found.and.fg_rank.eq.0)
741 c & write(iout,'(a15,f12.2,f8.1,2i5)')
742 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
750 if (newihpb(i).eq.idssb(j) .and.
751 & newjhpb(i).eq.jdssb(j)) found=.true.
755 c if (.not.found.and.fg_rank.eq.0)
756 c & write(iout,'(a15,f12.2,f8.1,2i5)')
757 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
772 c$$$c-----------------------------------------------------------------------------
774 c$$$ subroutine ss_relax(i_in,j_in)
778 c$$$ include 'DIMENSIONS'
779 c$$$ include 'COMMON.VAR'
780 c$$$ include 'COMMON.CHAIN'
781 c$$$ include 'COMMON.IOUNITS'
782 c$$$ include 'COMMON.INTERACT'
784 c$$$c Input arguments
785 c$$$ integer i_in,j_in
787 c$$$c Local variables
788 c$$$ integer i,iretcode,nfun_sc
790 c$$$ double precision var(maxvar),e_sc,etot
797 c$$$ mask_side(i_in)=1
798 c$$$ mask_side(j_in)=1
800 c$$$c Minimize the two selected side-chains
801 c$$$ call overlap_sc(scfail) ! Better not fail!
802 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
809 c$$$c-------------------------------------------------------------
811 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
812 c$$$c Minimize side-chains only, starting from geom but without modifying
814 c$$$c If mask_r is already set, only the selected side-chains are minimized,
815 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
819 c$$$ include 'DIMENSIONS'
820 c$$$ include 'COMMON.IOUNITS'
821 c$$$ include 'COMMON.VAR'
822 c$$$ include 'COMMON.CHAIN'
823 c$$$ include 'COMMON.GEO'
824 c$$$ include 'COMMON.MINIM'
826 c$$$ common /srutu/ icall
828 c$$$c Output arguments
829 c$$$ double precision etot_sc
830 c$$$ integer iretcode,nfun
832 c$$$c External functions/subroutines
833 c$$$ external func_sc,grad_sc,fdum
835 c$$$c Local variables
837 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
839 c$$$ double precision rdum(1)
840 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
842 c$$$ integer i,nvar_restr
845 c$$$cmc start_minim=.true.
846 c$$$ call deflt(2,iv,liv,lv,v)
847 c$$$* 12 means fresh start, dont call deflt
849 c$$$* max num of fun calls
850 c$$$ if (maxfun.eq.0) maxfun=500
852 c$$$* max num of iterations
853 c$$$ if (maxmin.eq.0) maxmin=1000
855 c$$$* controls output
857 c$$$* selects output unit
859 c$$$c iv(21)=iout ! DEBUG
860 c$$$c iv(21)=8 ! DEBUG
861 c$$$* 1 means to print out result
863 c$$$c iv(22)=1 ! DEBUG
864 c$$$* 1 means to print out summary stats
866 c$$$c iv(23)=1 ! DEBUG
867 c$$$* 1 means to print initial x and d
869 c$$$c iv(24)=1 ! DEBUG
870 c$$$* min val for v(radfac) default is 0.1
872 c$$$* max val for v(radfac) default is 4.0
875 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
876 c$$$* the sumsl default is 0.1
878 c$$$* false conv if (act fnctn decrease) .lt. v(34)
879 c$$$* the sumsl default is 100*machep
880 c$$$ v(34)=v(34)/100.0D0
881 c$$$* absolute convergence
882 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
884 c$$$* relative convergence
885 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
887 c$$$* controls initial step size
889 c$$$* large vals of d correspond to small components of step
893 c$$$ do i=nphi+1,nvar
897 c$$$ call geom_to_var(nvar,x)
898 c$$$ IF (mask_r) THEN
899 c$$$ do i=1,nres ! Just in case...
903 c$$$ call x2xx(x,xx,nvar_restr)
904 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
905 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
908 c$$$c When minimizing ALL side-chains, etotal_sc is a little
909 c$$$c faster if we don't set mask_r
915 c$$$ call x2xx(x,xx,nvar_restr)
916 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
917 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
920 c$$$ call var_to_geom(nvar,x)
921 c$$$ call chainbuild_sc
928 c$$$C--------------------------------------------------------------------------
930 c$$$ subroutine chainbuild_sc
932 c$$$ include 'DIMENSIONS'
933 c$$$ include 'COMMON.VAR'
934 c$$$ include 'COMMON.INTERACT'
936 c$$$c Local variables
941 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
942 c$$$ call locate_side_chain(i)
949 c$$$C--------------------------------------------------------------------------
951 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
955 c$$$ include 'DIMENSIONS'
956 c$$$ include 'COMMON.DERIV'
957 c$$$ include 'COMMON.VAR'
958 c$$$ include 'COMMON.MINIM'
959 c$$$ include 'COMMON.IOUNITS'
961 c$$$c Input arguments
963 c$$$ double precision x(maxvar)
964 c$$$ double precision ufparm
967 c$$$c Input/Output arguments
969 c$$$ integer uiparm(1)
970 c$$$ double precision urparm(1)
972 c$$$c Output arguments
973 c$$$ double precision f
975 c$$$c Local variables
976 c$$$ double precision energia(0:n_ene)
978 c$$$c Variables used to intercept NaNs
979 c$$$ double precision x_sum
988 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
991 c$$$ x_sum=x_sum+x(i_NAN)
993 c$$$c Calculate the energy only if the coordinates are ok
994 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
995 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
1001 c$$$ call var_to_geom_restr(n,x)
1003 c$$$ call chainbuild_sc
1004 c$$$ call etotal_sc(energia(0))
1006 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
1015 c$$$c-------------------------------------------------------
1017 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
1021 c$$$ include 'DIMENSIONS'
1022 c$$$ include 'COMMON.CHAIN'
1023 c$$$ include 'COMMON.DERIV'
1024 c$$$ include 'COMMON.VAR'
1025 c$$$ include 'COMMON.INTERACT'
1026 c$$$ include 'COMMON.MINIM'
1028 c$$$c Input arguments
1030 c$$$ double precision x(maxvar)
1031 c$$$ double precision ufparm
1032 c$$$ external ufparm
1034 c$$$c Input/Output arguments
1036 c$$$ integer uiparm(1)
1037 c$$$ double precision urparm(1)
1039 c$$$c Output arguments
1040 c$$$ double precision g(maxvar)
1042 c$$$c Local variables
1043 c$$$ double precision f,gphii,gthetai,galphai,gomegai
1044 c$$$ integer ig,ind,i,j,k,igall,ij
1047 c$$$ icg=mod(nf,2)+1
1048 c$$$ if (nf-nfl+1) 20,30,40
1049 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1050 c$$$c write (iout,*) 'grad 20'
1051 c$$$ if (nf.eq.0) return
1053 c$$$ 30 call var_to_geom_restr(n,x)
1054 c$$$ call chainbuild_sc
1056 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1058 c$$$ 40 call cartder
1060 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1066 c$$$ IF (mask_phi(i+2).eq.1) THEN
1068 c$$$ do j=i+1,nres-1
1071 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1072 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1078 c$$$ ind=ind+nres-1-i
1085 c$$$ IF (mask_theta(i+2).eq.1) THEN
1088 c$$$ do j=i+1,nres-1
1091 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1092 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1097 c$$$ ind=ind+nres-1-i
1102 c$$$ if (itype(i).ne.10) then
1103 c$$$ IF (mask_side(i).eq.1) THEN
1107 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1116 c$$$ if (itype(i).ne.10) then
1117 c$$$ IF (mask_side(i).eq.1) THEN
1121 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1129 c$$$C Add the components corresponding to local energy terms.
1136 c$$$ if (mask_phi(i).eq.1) then
1138 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1144 c$$$ if (mask_theta(i).eq.1) then
1146 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1152 c$$$ if (itype(i).ne.10) then
1154 c$$$ if (mask_side(i).eq.1) then
1156 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1163 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1169 c$$$C-----------------------------------------------------------------------------
1171 c$$$ subroutine etotal_sc(energy_sc)
1175 c$$$ include 'DIMENSIONS'
1176 c$$$ include 'COMMON.VAR'
1177 c$$$ include 'COMMON.INTERACT'
1178 c$$$ include 'COMMON.DERIV'
1179 c$$$ include 'COMMON.FFIELD'
1181 c$$$c Output arguments
1182 c$$$ double precision energy_sc(0:n_ene)
1184 c$$$c Local variables
1185 c$$$ double precision evdw,escloc
1190 c$$$ energy_sc(i)=0.0D0
1193 c$$$ if (mask_r) then
1194 c$$$ call egb_sc(evdw)
1195 c$$$ call esc_sc(escloc)
1198 c$$$ call esc(escloc)
1201 c$$$ if (evdw.eq.1.0D20) then
1202 c$$$ energy_sc(0)=evdw
1204 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1206 c$$$ energy_sc(1)=evdw
1207 c$$$ energy_sc(12)=escloc
1210 c$$$C Sum up the components of the Cartesian gradient.
1214 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1221 c$$$C-----------------------------------------------------------------------------
1223 c$$$ subroutine egb_sc(evdw)
1225 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1226 c$$$C assuming the Gay-Berne potential of interaction.
1228 c$$$ implicit real*8 (a-h,o-z)
1229 c$$$ include 'DIMENSIONS'
1230 c$$$ include 'COMMON.GEO'
1231 c$$$ include 'COMMON.VAR'
1232 c$$$ include 'COMMON.LOCAL'
1233 c$$$ include 'COMMON.CHAIN'
1234 c$$$ include 'COMMON.DERIV'
1235 c$$$ include 'COMMON.NAMES'
1236 c$$$ include 'COMMON.INTERACT'
1237 c$$$ include 'COMMON.IOUNITS'
1238 c$$$ include 'COMMON.CALC'
1239 c$$$ include 'COMMON.CONTROL'
1242 c$$$ energy_dec=.false.
1243 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1246 c$$$c if (icall.eq.0) lprn=.false.
1248 c$$$ do i=iatsc_s,iatsc_e
1250 c$$$ itypi1=itype(i+1)
1254 c$$$ dxi=dc_norm(1,nres+i)
1255 c$$$ dyi=dc_norm(2,nres+i)
1256 c$$$ dzi=dc_norm(3,nres+i)
1257 c$$$c dsci_inv=dsc_inv(itypi)
1258 c$$$ dsci_inv=vbld_inv(i+nres)
1259 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1260 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1262 c$$$C Calculate SC interaction energy.
1264 c$$$ do iint=1,nint_gr(i)
1265 c$$$ do j=istart(i,iint),iend(i,iint)
1266 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1269 c$$$c dscj_inv=dsc_inv(itypj)
1270 c$$$ dscj_inv=vbld_inv(j+nres)
1271 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1272 c$$$c & 1.0d0/vbld(j+nres)
1273 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1274 c$$$ sig0ij=sigma(itypi,itypj)
1275 c$$$ chi1=chi(itypi,itypj)
1276 c$$$ chi2=chi(itypj,itypi)
1277 c$$$ chi12=chi1*chi2
1278 c$$$ chip1=chip(itypi)
1279 c$$$ chip2=chip(itypj)
1280 c$$$ chip12=chip1*chip2
1281 c$$$ alf1=alp(itypi)
1282 c$$$ alf2=alp(itypj)
1283 c$$$ alf12=0.5D0*(alf1+alf2)
1284 c$$$C For diagnostics only!!!
1294 c$$$ xj=c(1,nres+j)-xi
1295 c$$$ yj=c(2,nres+j)-yi
1296 c$$$ zj=c(3,nres+j)-zi
1297 c$$$ dxj=dc_norm(1,nres+j)
1298 c$$$ dyj=dc_norm(2,nres+j)
1299 c$$$ dzj=dc_norm(3,nres+j)
1300 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1301 c$$$c write (iout,*) "j",j," dc_norm",
1302 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1303 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1304 c$$$ rij=dsqrt(rrij)
1305 c$$$C Calculate angle-dependent terms of energy and contributions to their
1307 c$$$ call sc_angular
1308 c$$$ sigsq=1.0D0/sigsq
1309 c$$$ sig=sig0ij*dsqrt(sigsq)
1310 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1311 c$$$c for diagnostics; uncomment
1312 c$$$c rij_shift=1.2*sig0ij
1313 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1314 c$$$ if (rij_shift.le.0.0D0) then
1316 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1317 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1318 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1321 c$$$ sigder=-sig*sigsq
1322 c$$$c---------------------------------------------------------------
1323 c$$$ rij_shift=1.0D0/rij_shift
1324 c$$$ fac=rij_shift**expon
1325 c$$$ e1=fac*fac*aa(itypi,itypj)
1326 c$$$ e2=fac*bb(itypi,itypj)
1327 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1328 c$$$ eps2der=evdwij*eps3rt
1329 c$$$ eps3der=evdwij*eps2rt
1330 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1331 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1332 c$$$ evdwij=evdwij*eps2rt*eps3rt
1333 c$$$ evdw=evdw+evdwij
1335 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1336 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1337 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1338 c$$$ & restyp(itypi),i,restyp(itypj),j,
1339 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1340 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1341 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1345 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1346 c$$$ & 'evdw',i,j,evdwij
1348 c$$$C Calculate gradient components.
1349 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1350 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1351 c$$$ sigder=fac*sigder
1354 c$$$C Calculate the radial part of the gradient
1358 c$$$C Calculate angular part of the gradient.
1364 c$$$ energy_dec=.false.
1368 c$$$c-----------------------------------------------------------------------------
1370 c$$$ subroutine esc_sc(escloc)
1371 c$$$C Calculate the local energy of a side chain and its derivatives in the
1372 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1373 c$$$C ALPHA and OMEGA.
1374 c$$$ implicit real*8 (a-h,o-z)
1375 c$$$ include 'DIMENSIONS'
1376 c$$$ include 'COMMON.GEO'
1377 c$$$ include 'COMMON.LOCAL'
1378 c$$$ include 'COMMON.VAR'
1379 c$$$ include 'COMMON.INTERACT'
1380 c$$$ include 'COMMON.DERIV'
1381 c$$$ include 'COMMON.CHAIN'
1382 c$$$ include 'COMMON.IOUNITS'
1383 c$$$ include 'COMMON.NAMES'
1384 c$$$ include 'COMMON.FFIELD'
1385 c$$$ include 'COMMON.CONTROL'
1386 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1387 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1388 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1389 c$$$ delta=0.02d0*pi
1391 c$$$c write (iout,'(a)') 'ESC'
1392 c$$$ do i=loc_start,loc_end
1393 c$$$ IF (mask_side(i).eq.1) THEN
1395 c$$$ if (it.eq.10) goto 1
1396 c$$$ nlobit=nlob(it)
1397 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1398 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1399 c$$$ theti=theta(i+1)-pipol
1400 c$$$ x(1)=dtan(theti)
1404 c$$$ if (x(2).gt.pi-delta) then
1406 c$$$ xtemp(2)=pi-delta
1408 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1410 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1411 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1412 c$$$ & escloci,dersc(2))
1413 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1414 c$$$ & ddersc0(1),dersc(1))
1415 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1416 c$$$ & ddersc0(3),dersc(3))
1417 c$$$ xtemp(2)=pi-delta
1418 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1420 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1421 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1422 c$$$ & dersc0(2),esclocbi,dersc02)
1423 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1424 c$$$ & dersc12,dersc01)
1425 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1426 c$$$ dersc0(1)=dersc01
1427 c$$$ dersc0(2)=dersc02
1428 c$$$ dersc0(3)=0.0d0
1430 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1432 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1433 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1434 c$$$c & esclocbi,ss,ssd
1435 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1436 c$$$c escloci=esclocbi
1437 c$$$c write (iout,*) escloci
1438 c$$$ else if (x(2).lt.delta) then
1442 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1444 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1445 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1446 c$$$ & escloci,dersc(2))
1447 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1448 c$$$ & ddersc0(1),dersc(1))
1449 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1450 c$$$ & ddersc0(3),dersc(3))
1452 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1454 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1455 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1456 c$$$ & dersc0(2),esclocbi,dersc02)
1457 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1458 c$$$ & dersc12,dersc01)
1459 c$$$ dersc0(1)=dersc01
1460 c$$$ dersc0(2)=dersc02
1461 c$$$ dersc0(3)=0.0d0
1462 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1464 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1466 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1467 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1468 c$$$c & esclocbi,ss,ssd
1469 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1470 c$$$c write (iout,*) escloci
1472 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1475 c$$$ escloc=escloc+escloci
1476 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1477 c$$$ & 'escloc',i,escloci
1478 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1480 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1481 c$$$ & wscloc*dersc(1)
1482 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1483 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1490 c$$$C-----------------------------------------------------------------------------
1492 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1494 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1495 c$$$C assuming the Gay-Berne potential of interaction.
1497 c$$$ implicit real*8 (a-h,o-z)
1498 c$$$ include 'DIMENSIONS'
1499 c$$$ include 'COMMON.GEO'
1500 c$$$ include 'COMMON.VAR'
1501 c$$$ include 'COMMON.LOCAL'
1502 c$$$ include 'COMMON.CHAIN'
1503 c$$$ include 'COMMON.DERIV'
1504 c$$$ include 'COMMON.NAMES'
1505 c$$$ include 'COMMON.INTERACT'
1506 c$$$ include 'COMMON.IOUNITS'
1507 c$$$ include 'COMMON.CALC'
1508 c$$$ include 'COMMON.CONTROL'
1511 c$$$ energy_dec=.false.
1512 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1516 c$$$c$$$ do i=iatsc_s,iatsc_e
1519 c$$$ itypi1=itype(i+1)
1523 c$$$ dxi=dc_norm(1,nres+i)
1524 c$$$ dyi=dc_norm(2,nres+i)
1525 c$$$ dzi=dc_norm(3,nres+i)
1526 c$$$c dsci_inv=dsc_inv(itypi)
1527 c$$$ dsci_inv=vbld_inv(i+nres)
1528 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1529 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1531 c$$$C Calculate SC interaction energy.
1533 c$$$c$$$ do iint=1,nint_gr(i)
1534 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1538 c$$$c dscj_inv=dsc_inv(itypj)
1539 c$$$ dscj_inv=vbld_inv(j+nres)
1540 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1541 c$$$c & 1.0d0/vbld(j+nres)
1542 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1543 c$$$ sig0ij=sigma(itypi,itypj)
1544 c$$$ chi1=chi(itypi,itypj)
1545 c$$$ chi2=chi(itypj,itypi)
1546 c$$$ chi12=chi1*chi2
1547 c$$$ chip1=chip(itypi)
1548 c$$$ chip2=chip(itypj)
1549 c$$$ chip12=chip1*chip2
1550 c$$$ alf1=alp(itypi)
1551 c$$$ alf2=alp(itypj)
1552 c$$$ alf12=0.5D0*(alf1+alf2)
1553 c$$$C For diagnostics only!!!
1563 c$$$ xj=c(1,nres+j)-xi
1564 c$$$ yj=c(2,nres+j)-yi
1565 c$$$ zj=c(3,nres+j)-zi
1566 c$$$ dxj=dc_norm(1,nres+j)
1567 c$$$ dyj=dc_norm(2,nres+j)
1568 c$$$ dzj=dc_norm(3,nres+j)
1569 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1570 c$$$c write (iout,*) "j",j," dc_norm",
1571 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1572 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1573 c$$$ rij=dsqrt(rrij)
1574 c$$$C Calculate angle-dependent terms of energy and contributions to their
1576 c$$$ call sc_angular
1577 c$$$ sigsq=1.0D0/sigsq
1578 c$$$ sig=sig0ij*dsqrt(sigsq)
1579 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1580 c$$$c for diagnostics; uncomment
1581 c$$$c rij_shift=1.2*sig0ij
1582 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1583 c$$$ if (rij_shift.le.0.0D0) then
1585 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1586 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1587 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1590 c$$$ sigder=-sig*sigsq
1591 c$$$c---------------------------------------------------------------
1592 c$$$ rij_shift=1.0D0/rij_shift
1593 c$$$ fac=rij_shift**expon
1594 c$$$ e1=fac*fac*aa(itypi,itypj)
1595 c$$$ e2=fac*bb(itypi,itypj)
1596 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1597 c$$$ eps2der=evdwij*eps3rt
1598 c$$$ eps3der=evdwij*eps2rt
1599 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1600 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1601 c$$$ evdwij=evdwij*eps2rt*eps3rt
1602 c$$$ evdw=evdw+evdwij
1604 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1605 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1606 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1607 c$$$ & restyp(itypi),i,restyp(itypj),j,
1608 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1609 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1610 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1614 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1615 c$$$ & 'evdw',i,j,evdwij
1617 c$$$C Calculate gradient components.
1618 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1619 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1620 c$$$ sigder=fac*sigder
1623 c$$$C Calculate the radial part of the gradient
1627 c$$$C Calculate angular part of the gradient.
1630 c$$$c$$$ enddo ! iint
1632 c$$$ energy_dec=.false.
1636 c$$$C-----------------------------------------------------------------------------
1638 c$$$ subroutine perturb_side_chain(i,angle)
1642 c$$$ include 'DIMENSIONS'
1643 c$$$ include 'COMMON.CHAIN'
1644 c$$$ include 'COMMON.GEO'
1645 c$$$ include 'COMMON.VAR'
1646 c$$$ include 'COMMON.LOCAL'
1647 c$$$ include 'COMMON.IOUNITS'
1649 c$$$c External functions
1650 c$$$ external ran_number
1651 c$$$ double precision ran_number
1653 c$$$c Input arguments
1655 c$$$ double precision angle ! In degrees
1657 c$$$c Local variables
1659 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1663 c$$$ rad_ang=angle*deg2rad
1666 c$$$ do while (length.lt.0.01)
1667 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1668 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1669 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1670 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1671 c$$$ + rand_v(3)*rand_v(3)
1672 c$$$ length=sqrt(length)
1673 c$$$ rand_v(1)=rand_v(1)/length
1674 c$$$ rand_v(2)=rand_v(2)/length
1675 c$$$ rand_v(3)=rand_v(3)/length
1676 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1677 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1678 c$$$ length=1.0D0-cost*cost
1679 c$$$ if (length.lt.0.0D0) length=0.0D0
1680 c$$$ length=sqrt(length)
1681 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1682 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1683 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1685 c$$$ rand_v(1)=rand_v(1)/length
1686 c$$$ rand_v(2)=rand_v(2)/length
1687 c$$$ rand_v(3)=rand_v(3)/length
1689 c$$$ cost=dcos(rad_ang)
1690 c$$$ sint=dsin(rad_ang)
1691 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1692 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1693 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1694 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1695 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1696 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1697 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1698 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1699 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1701 c$$$ call chainbuild_cart
1706 c$$$c----------------------------------------------------------------------------
1708 c$$$ subroutine ss_relax3(i_in,j_in)
1712 c$$$ include 'DIMENSIONS'
1713 c$$$ include 'COMMON.VAR'
1714 c$$$ include 'COMMON.CHAIN'
1715 c$$$ include 'COMMON.IOUNITS'
1716 c$$$ include 'COMMON.INTERACT'
1718 c$$$c External functions
1719 c$$$ external ran_number
1720 c$$$ double precision ran_number
1722 c$$$c Input arguments
1723 c$$$ integer i_in,j_in
1725 c$$$c Local variables
1726 c$$$ double precision energy_sc(0:n_ene),etot
1727 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1728 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1729 c$$$ integer n,i_pert,i
1730 c$$$ logical notdone
1739 c$$$ mask_side(i_in)=1
1740 c$$$ mask_side(j_in)=1
1742 c$$$ call etotal_sc(energy_sc)
1743 c$$$ etot=energy_sc(0)
1744 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1745 c$$$c + energy_sc(1),energy_sc(12)
1749 c$$$ do while (notdone)
1750 c$$$ if (mod(n,2).eq.0) then
1758 c$$$ org_dc(i)=dc(i,i_pert+nres)
1759 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1760 c$$$ org_c(i)=c(i,i_pert+nres)
1762 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1763 c$$$ call perturb_side_chain(i_pert,ang_pert)
1764 c$$$ call etotal_sc(energy_sc)
1765 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1766 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1767 c$$$ if (rand_fact.lt.exp_fact) then
1768 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1769 c$$$c + energy_sc(1),energy_sc(12)
1770 c$$$ etot=energy_sc(0)
1772 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1773 c$$$c + energy_sc(1),energy_sc(12)
1775 c$$$ dc(i,i_pert+nres)=org_dc(i)
1776 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1777 c$$$ c(i,i_pert+nres)=org_c(i)
1781 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1789 c$$$c----------------------------------------------------------------------------
1791 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1793 c$$$ include 'DIMENSIONS'
1795 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1796 c$$$*********************************************************************
1797 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1798 c$$$* the calling subprogram. *
1799 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1800 c$$$* calculated in the usual pythagorean way. *
1801 c$$$* absolute convergence occurs when the function is within v(31) of *
1802 c$$$* zero. unless you know the minimum value in advance, abs convg *
1803 c$$$* is probably not useful. *
1804 c$$$* relative convergence is when the model predicts that the function *
1805 c$$$* will decrease by less than v(32)*abs(fun). *
1806 c$$$*********************************************************************
1807 c$$$ include 'COMMON.IOUNITS'
1808 c$$$ include 'COMMON.VAR'
1809 c$$$ include 'COMMON.GEO'
1810 c$$$ include 'COMMON.MINIM'
1811 c$$$ include 'COMMON.CHAIN'
1813 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1814 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1815 c$$$ + orig_ss_dist(maxres2,maxres2)
1817 c$$$ double precision etot
1818 c$$$ integer iretcode,nfun,i_in,j_in
1821 c$$$ double precision dist
1822 c$$$ external ss_func,fdum
1823 c$$$ double precision ss_func,fdum
1825 c$$$ integer iv(liv),uiparm(2)
1826 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1830 c$$$ call deflt(2,iv,liv,lv,v)
1831 c$$$* 12 means fresh start, dont call deflt
1833 c$$$* max num of fun calls
1834 c$$$ if (maxfun.eq.0) maxfun=500
1836 c$$$* max num of iterations
1837 c$$$ if (maxmin.eq.0) maxmin=1000
1839 c$$$* controls output
1841 c$$$* selects output unit
1844 c$$$* 1 means to print out result
1846 c$$$* 1 means to print out summary stats
1848 c$$$* 1 means to print initial x and d
1850 c$$$* min val for v(radfac) default is 0.1
1852 c$$$* max val for v(radfac) default is 4.0
1855 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1856 c$$$* the sumsl default is 0.1
1858 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1859 c$$$* the sumsl default is 100*machep
1860 c$$$ v(34)=v(34)/100.0D0
1861 c$$$* absolute convergence
1862 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1865 c$$$* relative convergence
1866 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1869 c$$$* controls initial step size
1871 c$$$* large vals of d correspond to small components of step
1878 c$$$ orig_ss_dc(j,i)=dc(j,i)
1881 c$$$ call geom_to_var(nvar,orig_ss_var)
1885 c$$$ orig_ss_dist(j,i)=dist(j,i)
1886 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1887 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1888 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1900 c$$$ if (ialph(i,1).gt.0) then
1903 c$$$ x(k)=dc(j,i+nres)
1910 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1913 c$$$ nfun=iv(6)+iv(30)
1923 c$$$ if (ialph(i,1).gt.0) then
1926 c$$$ dc(j,i+nres)=x(k)
1930 c$$$ call chainbuild_cart
1935 c$$$C-----------------------------------------------------------------------------
1937 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1939 c$$$ include 'DIMENSIONS'
1940 c$$$ include 'COMMON.DERIV'
1941 c$$$ include 'COMMON.IOUNITS'
1942 c$$$ include 'COMMON.VAR'
1943 c$$$ include 'COMMON.CHAIN'
1944 c$$$ include 'COMMON.INTERACT'
1945 c$$$ include 'COMMON.SBRIDGE'
1947 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1948 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1949 c$$$ + orig_ss_dist(maxres2,maxres2)
1952 c$$$ double precision x(maxres6)
1954 c$$$ double precision f
1955 c$$$ integer uiparm(2)
1956 c$$$ real*8 urparm(1)
1957 c$$$ external ufparm
1958 c$$$ double precision ufparm
1961 c$$$ double precision dist
1963 c$$$ integer i,j,k,ss_i,ss_j
1964 c$$$ double precision tempf,var(maxvar)
1979 c$$$ if (ialph(i,1).gt.0) then
1982 c$$$ dc(j,i+nres)=x(k)
1986 c$$$ call chainbuild_cart
1988 c$$$ call geom_to_var(nvar,var)
1990 c$$$c Constraints on all angles
1992 c$$$ tempf=var(i)-orig_ss_var(i)
1993 c$$$ f=f+tempf*tempf
1996 c$$$c Constraints on all distances
1998 c$$$ if (i.gt.1) then
1999 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
2000 c$$$ f=f+tempf*tempf
2003 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
2004 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
2005 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
2006 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2007 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
2008 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2009 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
2010 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2014 c$$$c Constraints for the relevant CYS-CYS
2015 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
2016 c$$$ f=f+tempf*tempf
2017 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
2019 c$$$c$$$ if (nf.ne.nfl) then
2020 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
2021 c$$$c$$$ + f,dist(5+nres,14+nres)
2029 c$$$C-----------------------------------------------------------------------------
2030 c$$$C-----------------------------------------------------------------------------
2031 subroutine triple_ssbond_ene(resi,resj,resk,eij)
2032 implicit real*8 (a-h,o-z)
2033 include 'DIMENSIONS'
2034 include 'COMMON.SBRIDGE'
2035 include 'COMMON.CHAIN'
2036 include 'COMMON.DERIV'
2037 include 'COMMON.LOCAL'
2038 include 'COMMON.INTERACT'
2039 include 'COMMON.VAR'
2040 include 'COMMON.IOUNITS'
2041 include 'COMMON.CALC'
2044 C include 'COMMON.MD'
2048 c External functions
2049 double precision h_base
2053 integer resi,resj,resk
2056 double precision eij,eij1,eij2,eij3
2060 c integer itypi,itypj,k,l
2061 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2062 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2063 double precision xik,yik,zik,xjk,yjk,zjk
2064 double precision sig0ij,ljd,sig,fac,e1,e2
2065 double precision dcosom1(3),dcosom2(3),ed
2066 double precision pom1,pom2
2067 double precision ljA,ljB,ljXs
2068 double precision d_ljB(1:3)
2069 double precision ssA,ssB,ssC,ssXs
2070 double precision ssxm,ljxm,ssm,ljm
2071 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2076 C write(iout,*) resi,resj,resk
2078 dxi=dc_norm(1,nres+i)
2079 dyi=dc_norm(2,nres+i)
2080 dzi=dc_norm(3,nres+i)
2081 dsci_inv=vbld_inv(i+nres)
2091 dxj=dc_norm(1,nres+j)
2092 dyj=dc_norm(2,nres+j)
2093 dzj=dc_norm(3,nres+j)
2094 dscj_inv=vbld_inv(j+nres)
2100 dxk=dc_norm(1,nres+k)
2101 dyk=dc_norm(2,nres+k)
2102 dzk=dc_norm(3,nres+k)
2103 dscj_inv=vbld_inv(k+nres)
2113 rrij=(xij*xij+yij*yij+zij*zij)
2114 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2115 rrik=(xik*xik+yik*yik+zik*zik)
2117 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2119 C there are three combination of distances for each trisulfide bonds
2120 C The first case the ith atom is the center
2121 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2122 C distance y is second distance the a,b,c,d are parameters derived for
2123 C this problem d parameter was set as a penalty currenlty set to 1.
2124 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2125 C second case jth atom is center
2126 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2127 C the third case kth atom is the center
2128 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2133 C write(iout,*)i,j,k,eij
2134 C The energy penalty calculated now time for the gradient part
2135 C derivative over rij
2136 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2137 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2142 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2143 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2146 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2147 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2149 C now derivative over rik
2150 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2151 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2156 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2157 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2160 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2161 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2163 C now derivative over rjk
2164 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2165 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2170 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2171 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2174 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2175 gvdwc(l,k)=gvdwc(l,k)+gg(l)