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----------------------------------------------------------------------------
730 subroutine read_ssHist
735 include "DIMENSIONS.FREE"
736 include 'COMMON.FREE'
740 character*80 controlcard
743 call card_concat(controlcard,.true.)
745 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
752 c----------------------------------------------------------------------------
755 C-----------------------------------------------------------------------------
756 C-----------------------------------------------------------------------------
757 C-----------------------------------------------------------------------------
758 C-----------------------------------------------------------------------------
759 C-----------------------------------------------------------------------------
760 C-----------------------------------------------------------------------------
761 C-----------------------------------------------------------------------------
763 c$$$c-----------------------------------------------------------------------------
765 c$$$ subroutine ss_relax(i_in,j_in)
769 c$$$ include 'DIMENSIONS'
770 c$$$ include 'COMMON.VAR'
771 c$$$ include 'COMMON.CHAIN'
772 c$$$ include 'COMMON.IOUNITS'
773 c$$$ include 'COMMON.INTERACT'
775 c$$$c Input arguments
776 c$$$ integer i_in,j_in
778 c$$$c Local variables
779 c$$$ integer i,iretcode,nfun_sc
781 c$$$ double precision var(maxvar),e_sc,etot
788 c$$$ mask_side(i_in)=1
789 c$$$ mask_side(j_in)=1
791 c$$$c Minimize the two selected side-chains
792 c$$$ call overlap_sc(scfail) ! Better not fail!
793 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
800 c$$$c-------------------------------------------------------------
802 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
803 c$$$c Minimize side-chains only, starting from geom but without modifying
805 c$$$c If mask_r is already set, only the selected side-chains are minimized,
806 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
810 c$$$ include 'DIMENSIONS'
811 c$$$ include 'COMMON.IOUNITS'
812 c$$$ include 'COMMON.VAR'
813 c$$$ include 'COMMON.CHAIN'
814 c$$$ include 'COMMON.GEO'
815 c$$$ include 'COMMON.MINIM'
817 c$$$ common /srutu/ icall
819 c$$$c Output arguments
820 c$$$ double precision etot_sc
821 c$$$ integer iretcode,nfun
823 c$$$c External functions/subroutines
824 c$$$ external func_sc,grad_sc,fdum
826 c$$$c Local variables
828 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
830 c$$$ double precision rdum(1)
831 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
833 c$$$ integer i,nvar_restr
836 c$$$cmc start_minim=.true.
837 c$$$ call deflt(2,iv,liv,lv,v)
838 c$$$* 12 means fresh start, dont call deflt
840 c$$$* max num of fun calls
841 c$$$ if (maxfun.eq.0) maxfun=500
843 c$$$* max num of iterations
844 c$$$ if (maxmin.eq.0) maxmin=1000
846 c$$$* controls output
848 c$$$* selects output unit
850 c$$$c iv(21)=iout ! DEBUG
851 c$$$c iv(21)=8 ! DEBUG
852 c$$$* 1 means to print out result
854 c$$$c iv(22)=1 ! DEBUG
855 c$$$* 1 means to print out summary stats
857 c$$$c iv(23)=1 ! DEBUG
858 c$$$* 1 means to print initial x and d
860 c$$$c iv(24)=1 ! DEBUG
861 c$$$* min val for v(radfac) default is 0.1
863 c$$$* max val for v(radfac) default is 4.0
866 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
867 c$$$* the sumsl default is 0.1
869 c$$$* false conv if (act fnctn decrease) .lt. v(34)
870 c$$$* the sumsl default is 100*machep
871 c$$$ v(34)=v(34)/100.0D0
872 c$$$* absolute convergence
873 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
875 c$$$* relative convergence
876 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
878 c$$$* controls initial step size
880 c$$$* large vals of d correspond to small components of step
884 c$$$ do i=nphi+1,nvar
888 c$$$ call geom_to_var(nvar,x)
889 c$$$ IF (mask_r) THEN
890 c$$$ do i=1,nres ! Just in case...
894 c$$$ call x2xx(x,xx,nvar_restr)
895 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
896 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
899 c$$$c When minimizing ALL side-chains, etotal_sc is a little
900 c$$$c faster if we don't set mask_r
906 c$$$ call x2xx(x,xx,nvar_restr)
907 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
908 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
911 c$$$ call var_to_geom(nvar,x)
912 c$$$ call chainbuild_sc
919 c$$$C--------------------------------------------------------------------------
921 c$$$ subroutine chainbuild_sc
923 c$$$ include 'DIMENSIONS'
924 c$$$ include 'COMMON.VAR'
925 c$$$ include 'COMMON.INTERACT'
927 c$$$c Local variables
932 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
933 c$$$ call locate_side_chain(i)
940 c$$$C--------------------------------------------------------------------------
942 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
946 c$$$ include 'DIMENSIONS'
947 c$$$ include 'COMMON.DERIV'
948 c$$$ include 'COMMON.VAR'
949 c$$$ include 'COMMON.MINIM'
950 c$$$ include 'COMMON.IOUNITS'
952 c$$$c Input arguments
954 c$$$ double precision x(maxvar)
955 c$$$ double precision ufparm
958 c$$$c Input/Output arguments
960 c$$$ integer uiparm(1)
961 c$$$ double precision urparm(1)
963 c$$$c Output arguments
964 c$$$ double precision f
966 c$$$c Local variables
967 c$$$ double precision energia(0:n_ene)
969 c$$$c Variables used to intercept NaNs
970 c$$$ double precision x_sum
979 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
982 c$$$ x_sum=x_sum+x(i_NAN)
984 c$$$c Calculate the energy only if the coordinates are ok
985 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
986 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
992 c$$$ call var_to_geom_restr(n,x)
994 c$$$ call chainbuild_sc
995 c$$$ call etotal_sc(energia(0))
997 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
1006 c$$$c-------------------------------------------------------
1008 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
1012 c$$$ include 'DIMENSIONS'
1013 c$$$ include 'COMMON.CHAIN'
1014 c$$$ include 'COMMON.DERIV'
1015 c$$$ include 'COMMON.VAR'
1016 c$$$ include 'COMMON.INTERACT'
1017 c$$$ include 'COMMON.MINIM'
1019 c$$$c Input arguments
1021 c$$$ double precision x(maxvar)
1022 c$$$ double precision ufparm
1023 c$$$ external ufparm
1025 c$$$c Input/Output arguments
1027 c$$$ integer uiparm(1)
1028 c$$$ double precision urparm(1)
1030 c$$$c Output arguments
1031 c$$$ double precision g(maxvar)
1033 c$$$c Local variables
1034 c$$$ double precision f,gphii,gthetai,galphai,gomegai
1035 c$$$ integer ig,ind,i,j,k,igall,ij
1038 c$$$ icg=mod(nf,2)+1
1039 c$$$ if (nf-nfl+1) 20,30,40
1040 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1041 c$$$c write (iout,*) 'grad 20'
1042 c$$$ if (nf.eq.0) return
1044 c$$$ 30 call var_to_geom_restr(n,x)
1045 c$$$ call chainbuild_sc
1047 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1049 c$$$ 40 call cartder
1051 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1057 c$$$ IF (mask_phi(i+2).eq.1) THEN
1059 c$$$ do j=i+1,nres-1
1062 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1063 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1069 c$$$ ind=ind+nres-1-i
1076 c$$$ IF (mask_theta(i+2).eq.1) THEN
1079 c$$$ do j=i+1,nres-1
1082 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1083 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1088 c$$$ ind=ind+nres-1-i
1093 c$$$ if (itype(i).ne.10) then
1094 c$$$ IF (mask_side(i).eq.1) THEN
1098 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1107 c$$$ if (itype(i).ne.10) then
1108 c$$$ IF (mask_side(i).eq.1) THEN
1112 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1120 c$$$C Add the components corresponding to local energy terms.
1127 c$$$ if (mask_phi(i).eq.1) then
1129 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1135 c$$$ if (mask_theta(i).eq.1) then
1137 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1143 c$$$ if (itype(i).ne.10) then
1145 c$$$ if (mask_side(i).eq.1) then
1147 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1154 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1160 c$$$C-----------------------------------------------------------------------------
1162 c$$$ subroutine etotal_sc(energy_sc)
1166 c$$$ include 'DIMENSIONS'
1167 c$$$ include 'COMMON.VAR'
1168 c$$$ include 'COMMON.INTERACT'
1169 c$$$ include 'COMMON.DERIV'
1170 c$$$ include 'COMMON.FFIELD'
1172 c$$$c Output arguments
1173 c$$$ double precision energy_sc(0:n_ene)
1175 c$$$c Local variables
1176 c$$$ double precision evdw,escloc
1181 c$$$ energy_sc(i)=0.0D0
1184 c$$$ if (mask_r) then
1185 c$$$ call egb_sc(evdw)
1186 c$$$ call esc_sc(escloc)
1189 c$$$ call esc(escloc)
1192 c$$$ if (evdw.eq.1.0D20) then
1193 c$$$ energy_sc(0)=evdw
1195 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1197 c$$$ energy_sc(1)=evdw
1198 c$$$ energy_sc(12)=escloc
1201 c$$$C Sum up the components of the Cartesian gradient.
1205 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1212 c$$$C-----------------------------------------------------------------------------
1214 c$$$ subroutine egb_sc(evdw)
1216 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1217 c$$$C assuming the Gay-Berne potential of interaction.
1219 c$$$ implicit real*8 (a-h,o-z)
1220 c$$$ include 'DIMENSIONS'
1221 c$$$ include 'COMMON.GEO'
1222 c$$$ include 'COMMON.VAR'
1223 c$$$ include 'COMMON.LOCAL'
1224 c$$$ include 'COMMON.CHAIN'
1225 c$$$ include 'COMMON.DERIV'
1226 c$$$ include 'COMMON.NAMES'
1227 c$$$ include 'COMMON.INTERACT'
1228 c$$$ include 'COMMON.IOUNITS'
1229 c$$$ include 'COMMON.CALC'
1230 c$$$ include 'COMMON.CONTROL'
1233 c$$$ energy_dec=.false.
1234 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1237 c$$$c if (icall.eq.0) lprn=.false.
1239 c$$$ do i=iatsc_s,iatsc_e
1241 c$$$ itypi1=itype(i+1)
1245 c$$$ dxi=dc_norm(1,nres+i)
1246 c$$$ dyi=dc_norm(2,nres+i)
1247 c$$$ dzi=dc_norm(3,nres+i)
1248 c$$$c dsci_inv=dsc_inv(itypi)
1249 c$$$ dsci_inv=vbld_inv(i+nres)
1250 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1251 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1253 c$$$C Calculate SC interaction energy.
1255 c$$$ do iint=1,nint_gr(i)
1256 c$$$ do j=istart(i,iint),iend(i,iint)
1257 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1260 c$$$c dscj_inv=dsc_inv(itypj)
1261 c$$$ dscj_inv=vbld_inv(j+nres)
1262 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1263 c$$$c & 1.0d0/vbld(j+nres)
1264 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1265 c$$$ sig0ij=sigma(itypi,itypj)
1266 c$$$ chi1=chi(itypi,itypj)
1267 c$$$ chi2=chi(itypj,itypi)
1268 c$$$ chi12=chi1*chi2
1269 c$$$ chip1=chip(itypi)
1270 c$$$ chip2=chip(itypj)
1271 c$$$ chip12=chip1*chip2
1272 c$$$ alf1=alp(itypi)
1273 c$$$ alf2=alp(itypj)
1274 c$$$ alf12=0.5D0*(alf1+alf2)
1275 c$$$C For diagnostics only!!!
1285 c$$$ xj=c(1,nres+j)-xi
1286 c$$$ yj=c(2,nres+j)-yi
1287 c$$$ zj=c(3,nres+j)-zi
1288 c$$$ dxj=dc_norm(1,nres+j)
1289 c$$$ dyj=dc_norm(2,nres+j)
1290 c$$$ dzj=dc_norm(3,nres+j)
1291 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1292 c$$$c write (iout,*) "j",j," dc_norm",
1293 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1294 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1295 c$$$ rij=dsqrt(rrij)
1296 c$$$C Calculate angle-dependent terms of energy and contributions to their
1298 c$$$ call sc_angular
1299 c$$$ sigsq=1.0D0/sigsq
1300 c$$$ sig=sig0ij*dsqrt(sigsq)
1301 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1302 c$$$c for diagnostics; uncomment
1303 c$$$c rij_shift=1.2*sig0ij
1304 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1305 c$$$ if (rij_shift.le.0.0D0) then
1307 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1308 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1309 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1312 c$$$ sigder=-sig*sigsq
1313 c$$$c---------------------------------------------------------------
1314 c$$$ rij_shift=1.0D0/rij_shift
1315 c$$$ fac=rij_shift**expon
1316 c$$$ e1=fac*fac*aa(itypi,itypj)
1317 c$$$ e2=fac*bb(itypi,itypj)
1318 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1319 c$$$ eps2der=evdwij*eps3rt
1320 c$$$ eps3der=evdwij*eps2rt
1321 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1322 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1323 c$$$ evdwij=evdwij*eps2rt*eps3rt
1324 c$$$ evdw=evdw+evdwij
1326 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1327 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1328 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1329 c$$$ & restyp(itypi),i,restyp(itypj),j,
1330 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1331 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1332 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1336 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1337 c$$$ & 'evdw',i,j,evdwij
1339 c$$$C Calculate gradient components.
1340 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1341 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1342 c$$$ sigder=fac*sigder
1345 c$$$C Calculate the radial part of the gradient
1349 c$$$C Calculate angular part of the gradient.
1355 c$$$ energy_dec=.false.
1359 c$$$c-----------------------------------------------------------------------------
1361 c$$$ subroutine esc_sc(escloc)
1362 c$$$C Calculate the local energy of a side chain and its derivatives in the
1363 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1364 c$$$C ALPHA and OMEGA.
1365 c$$$ implicit real*8 (a-h,o-z)
1366 c$$$ include 'DIMENSIONS'
1367 c$$$ include 'COMMON.GEO'
1368 c$$$ include 'COMMON.LOCAL'
1369 c$$$ include 'COMMON.VAR'
1370 c$$$ include 'COMMON.INTERACT'
1371 c$$$ include 'COMMON.DERIV'
1372 c$$$ include 'COMMON.CHAIN'
1373 c$$$ include 'COMMON.IOUNITS'
1374 c$$$ include 'COMMON.NAMES'
1375 c$$$ include 'COMMON.FFIELD'
1376 c$$$ include 'COMMON.CONTROL'
1377 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1378 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1379 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1380 c$$$ delta=0.02d0*pi
1382 c$$$c write (iout,'(a)') 'ESC'
1383 c$$$ do i=loc_start,loc_end
1384 c$$$ IF (mask_side(i).eq.1) THEN
1386 c$$$ if (it.eq.10) goto 1
1387 c$$$ nlobit=nlob(it)
1388 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1389 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1390 c$$$ theti=theta(i+1)-pipol
1391 c$$$ x(1)=dtan(theti)
1395 c$$$ if (x(2).gt.pi-delta) then
1397 c$$$ xtemp(2)=pi-delta
1399 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1401 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1402 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1403 c$$$ & escloci,dersc(2))
1404 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1405 c$$$ & ddersc0(1),dersc(1))
1406 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1407 c$$$ & ddersc0(3),dersc(3))
1408 c$$$ xtemp(2)=pi-delta
1409 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1411 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1412 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1413 c$$$ & dersc0(2),esclocbi,dersc02)
1414 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1415 c$$$ & dersc12,dersc01)
1416 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1417 c$$$ dersc0(1)=dersc01
1418 c$$$ dersc0(2)=dersc02
1419 c$$$ dersc0(3)=0.0d0
1421 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1423 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1424 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1425 c$$$c & esclocbi,ss,ssd
1426 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1427 c$$$c escloci=esclocbi
1428 c$$$c write (iout,*) escloci
1429 c$$$ else if (x(2).lt.delta) then
1433 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1435 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1436 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1437 c$$$ & escloci,dersc(2))
1438 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1439 c$$$ & ddersc0(1),dersc(1))
1440 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1441 c$$$ & ddersc0(3),dersc(3))
1443 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1445 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1446 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1447 c$$$ & dersc0(2),esclocbi,dersc02)
1448 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1449 c$$$ & dersc12,dersc01)
1450 c$$$ dersc0(1)=dersc01
1451 c$$$ dersc0(2)=dersc02
1452 c$$$ dersc0(3)=0.0d0
1453 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1455 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1457 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1458 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1459 c$$$c & esclocbi,ss,ssd
1460 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1461 c$$$c write (iout,*) escloci
1463 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1466 c$$$ escloc=escloc+escloci
1467 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1468 c$$$ & 'escloc',i,escloci
1469 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1471 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1472 c$$$ & wscloc*dersc(1)
1473 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1474 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1481 c$$$C-----------------------------------------------------------------------------
1483 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1485 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1486 c$$$C assuming the Gay-Berne potential of interaction.
1488 c$$$ implicit real*8 (a-h,o-z)
1489 c$$$ include 'DIMENSIONS'
1490 c$$$ include 'COMMON.GEO'
1491 c$$$ include 'COMMON.VAR'
1492 c$$$ include 'COMMON.LOCAL'
1493 c$$$ include 'COMMON.CHAIN'
1494 c$$$ include 'COMMON.DERIV'
1495 c$$$ include 'COMMON.NAMES'
1496 c$$$ include 'COMMON.INTERACT'
1497 c$$$ include 'COMMON.IOUNITS'
1498 c$$$ include 'COMMON.CALC'
1499 c$$$ include 'COMMON.CONTROL'
1502 c$$$ energy_dec=.false.
1503 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1507 c$$$c$$$ do i=iatsc_s,iatsc_e
1510 c$$$ itypi1=itype(i+1)
1514 c$$$ dxi=dc_norm(1,nres+i)
1515 c$$$ dyi=dc_norm(2,nres+i)
1516 c$$$ dzi=dc_norm(3,nres+i)
1517 c$$$c dsci_inv=dsc_inv(itypi)
1518 c$$$ dsci_inv=vbld_inv(i+nres)
1519 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1520 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1522 c$$$C Calculate SC interaction energy.
1524 c$$$c$$$ do iint=1,nint_gr(i)
1525 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1529 c$$$c dscj_inv=dsc_inv(itypj)
1530 c$$$ dscj_inv=vbld_inv(j+nres)
1531 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1532 c$$$c & 1.0d0/vbld(j+nres)
1533 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1534 c$$$ sig0ij=sigma(itypi,itypj)
1535 c$$$ chi1=chi(itypi,itypj)
1536 c$$$ chi2=chi(itypj,itypi)
1537 c$$$ chi12=chi1*chi2
1538 c$$$ chip1=chip(itypi)
1539 c$$$ chip2=chip(itypj)
1540 c$$$ chip12=chip1*chip2
1541 c$$$ alf1=alp(itypi)
1542 c$$$ alf2=alp(itypj)
1543 c$$$ alf12=0.5D0*(alf1+alf2)
1544 c$$$C For diagnostics only!!!
1554 c$$$ xj=c(1,nres+j)-xi
1555 c$$$ yj=c(2,nres+j)-yi
1556 c$$$ zj=c(3,nres+j)-zi
1557 c$$$ dxj=dc_norm(1,nres+j)
1558 c$$$ dyj=dc_norm(2,nres+j)
1559 c$$$ dzj=dc_norm(3,nres+j)
1560 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1561 c$$$c write (iout,*) "j",j," dc_norm",
1562 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1563 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1564 c$$$ rij=dsqrt(rrij)
1565 c$$$C Calculate angle-dependent terms of energy and contributions to their
1567 c$$$ call sc_angular
1568 c$$$ sigsq=1.0D0/sigsq
1569 c$$$ sig=sig0ij*dsqrt(sigsq)
1570 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1571 c$$$c for diagnostics; uncomment
1572 c$$$c rij_shift=1.2*sig0ij
1573 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1574 c$$$ if (rij_shift.le.0.0D0) then
1576 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1577 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1578 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1581 c$$$ sigder=-sig*sigsq
1582 c$$$c---------------------------------------------------------------
1583 c$$$ rij_shift=1.0D0/rij_shift
1584 c$$$ fac=rij_shift**expon
1585 c$$$ e1=fac*fac*aa(itypi,itypj)
1586 c$$$ e2=fac*bb(itypi,itypj)
1587 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1588 c$$$ eps2der=evdwij*eps3rt
1589 c$$$ eps3der=evdwij*eps2rt
1590 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1591 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1592 c$$$ evdwij=evdwij*eps2rt*eps3rt
1593 c$$$ evdw=evdw+evdwij
1595 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1596 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1597 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1598 c$$$ & restyp(itypi),i,restyp(itypj),j,
1599 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1600 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1601 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1605 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1606 c$$$ & 'evdw',i,j,evdwij
1608 c$$$C Calculate gradient components.
1609 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1610 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1611 c$$$ sigder=fac*sigder
1614 c$$$C Calculate the radial part of the gradient
1618 c$$$C Calculate angular part of the gradient.
1621 c$$$c$$$ enddo ! iint
1623 c$$$ energy_dec=.false.
1627 c$$$C-----------------------------------------------------------------------------
1629 c$$$ subroutine perturb_side_chain(i,angle)
1633 c$$$ include 'DIMENSIONS'
1634 c$$$ include 'COMMON.CHAIN'
1635 c$$$ include 'COMMON.GEO'
1636 c$$$ include 'COMMON.VAR'
1637 c$$$ include 'COMMON.LOCAL'
1638 c$$$ include 'COMMON.IOUNITS'
1640 c$$$c External functions
1641 c$$$ external ran_number
1642 c$$$ double precision ran_number
1644 c$$$c Input arguments
1646 c$$$ double precision angle ! In degrees
1648 c$$$c Local variables
1650 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1654 c$$$ rad_ang=angle*deg2rad
1657 c$$$ do while (length.lt.0.01)
1658 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1659 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1660 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1661 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1662 c$$$ + rand_v(3)*rand_v(3)
1663 c$$$ length=sqrt(length)
1664 c$$$ rand_v(1)=rand_v(1)/length
1665 c$$$ rand_v(2)=rand_v(2)/length
1666 c$$$ rand_v(3)=rand_v(3)/length
1667 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1668 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1669 c$$$ length=1.0D0-cost*cost
1670 c$$$ if (length.lt.0.0D0) length=0.0D0
1671 c$$$ length=sqrt(length)
1672 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1673 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1674 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1676 c$$$ rand_v(1)=rand_v(1)/length
1677 c$$$ rand_v(2)=rand_v(2)/length
1678 c$$$ rand_v(3)=rand_v(3)/length
1680 c$$$ cost=dcos(rad_ang)
1681 c$$$ sint=dsin(rad_ang)
1682 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1683 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1684 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1685 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1686 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1687 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1688 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1689 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1690 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1692 c$$$ call chainbuild_cart
1697 c$$$c----------------------------------------------------------------------------
1699 c$$$ subroutine ss_relax3(i_in,j_in)
1703 c$$$ include 'DIMENSIONS'
1704 c$$$ include 'COMMON.VAR'
1705 c$$$ include 'COMMON.CHAIN'
1706 c$$$ include 'COMMON.IOUNITS'
1707 c$$$ include 'COMMON.INTERACT'
1709 c$$$c External functions
1710 c$$$ external ran_number
1711 c$$$ double precision ran_number
1713 c$$$c Input arguments
1714 c$$$ integer i_in,j_in
1716 c$$$c Local variables
1717 c$$$ double precision energy_sc(0:n_ene),etot
1718 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1719 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1720 c$$$ integer n,i_pert,i
1721 c$$$ logical notdone
1730 c$$$ mask_side(i_in)=1
1731 c$$$ mask_side(j_in)=1
1733 c$$$ call etotal_sc(energy_sc)
1734 c$$$ etot=energy_sc(0)
1735 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1736 c$$$c + energy_sc(1),energy_sc(12)
1740 c$$$ do while (notdone)
1741 c$$$ if (mod(n,2).eq.0) then
1749 c$$$ org_dc(i)=dc(i,i_pert+nres)
1750 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1751 c$$$ org_c(i)=c(i,i_pert+nres)
1753 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1754 c$$$ call perturb_side_chain(i_pert,ang_pert)
1755 c$$$ call etotal_sc(energy_sc)
1756 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1757 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1758 c$$$ if (rand_fact.lt.exp_fact) then
1759 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1760 c$$$c + energy_sc(1),energy_sc(12)
1761 c$$$ etot=energy_sc(0)
1763 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1764 c$$$c + energy_sc(1),energy_sc(12)
1766 c$$$ dc(i,i_pert+nres)=org_dc(i)
1767 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1768 c$$$ c(i,i_pert+nres)=org_c(i)
1772 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1780 c$$$c----------------------------------------------------------------------------
1782 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1784 c$$$ include 'DIMENSIONS'
1786 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1787 c$$$*********************************************************************
1788 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1789 c$$$* the calling subprogram. *
1790 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1791 c$$$* calculated in the usual pythagorean way. *
1792 c$$$* absolute convergence occurs when the function is within v(31) of *
1793 c$$$* zero. unless you know the minimum value in advance, abs convg *
1794 c$$$* is probably not useful. *
1795 c$$$* relative convergence is when the model predicts that the function *
1796 c$$$* will decrease by less than v(32)*abs(fun). *
1797 c$$$*********************************************************************
1798 c$$$ include 'COMMON.IOUNITS'
1799 c$$$ include 'COMMON.VAR'
1800 c$$$ include 'COMMON.GEO'
1801 c$$$ include 'COMMON.MINIM'
1802 c$$$ include 'COMMON.CHAIN'
1804 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1805 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1806 c$$$ + orig_ss_dist(maxres2,maxres2)
1808 c$$$ double precision etot
1809 c$$$ integer iretcode,nfun,i_in,j_in
1812 c$$$ double precision dist
1813 c$$$ external ss_func,fdum
1814 c$$$ double precision ss_func,fdum
1816 c$$$ integer iv(liv),uiparm(2)
1817 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1821 c$$$ call deflt(2,iv,liv,lv,v)
1822 c$$$* 12 means fresh start, dont call deflt
1824 c$$$* max num of fun calls
1825 c$$$ if (maxfun.eq.0) maxfun=500
1827 c$$$* max num of iterations
1828 c$$$ if (maxmin.eq.0) maxmin=1000
1830 c$$$* controls output
1832 c$$$* selects output unit
1835 c$$$* 1 means to print out result
1837 c$$$* 1 means to print out summary stats
1839 c$$$* 1 means to print initial x and d
1841 c$$$* min val for v(radfac) default is 0.1
1843 c$$$* max val for v(radfac) default is 4.0
1846 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1847 c$$$* the sumsl default is 0.1
1849 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1850 c$$$* the sumsl default is 100*machep
1851 c$$$ v(34)=v(34)/100.0D0
1852 c$$$* absolute convergence
1853 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1856 c$$$* relative convergence
1857 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1860 c$$$* controls initial step size
1862 c$$$* large vals of d correspond to small components of step
1869 c$$$ orig_ss_dc(j,i)=dc(j,i)
1872 c$$$ call geom_to_var(nvar,orig_ss_var)
1876 c$$$ orig_ss_dist(j,i)=dist(j,i)
1877 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1878 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1879 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1891 c$$$ if (ialph(i,1).gt.0) then
1894 c$$$ x(k)=dc(j,i+nres)
1901 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1904 c$$$ nfun=iv(6)+iv(30)
1914 c$$$ if (ialph(i,1).gt.0) then
1917 c$$$ dc(j,i+nres)=x(k)
1921 c$$$ call chainbuild_cart
1926 c$$$C-----------------------------------------------------------------------------
1928 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1930 c$$$ include 'DIMENSIONS'
1931 c$$$ include 'COMMON.DERIV'
1932 c$$$ include 'COMMON.IOUNITS'
1933 c$$$ include 'COMMON.VAR'
1934 c$$$ include 'COMMON.CHAIN'
1935 c$$$ include 'COMMON.INTERACT'
1936 c$$$ include 'COMMON.SBRIDGE'
1938 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1939 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1940 c$$$ + orig_ss_dist(maxres2,maxres2)
1943 c$$$ double precision x(maxres6)
1945 c$$$ double precision f
1946 c$$$ integer uiparm(2)
1947 c$$$ real*8 urparm(1)
1948 c$$$ external ufparm
1949 c$$$ double precision ufparm
1952 c$$$ double precision dist
1954 c$$$ integer i,j,k,ss_i,ss_j
1955 c$$$ double precision tempf,var(maxvar)
1970 c$$$ if (ialph(i,1).gt.0) then
1973 c$$$ dc(j,i+nres)=x(k)
1977 c$$$ call chainbuild_cart
1979 c$$$ call geom_to_var(nvar,var)
1981 c$$$c Constraints on all angles
1983 c$$$ tempf=var(i)-orig_ss_var(i)
1984 c$$$ f=f+tempf*tempf
1987 c$$$c Constraints on all distances
1989 c$$$ if (i.gt.1) then
1990 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1991 c$$$ f=f+tempf*tempf
1994 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1995 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1996 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1997 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1998 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1999 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2000 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
2001 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2005 c$$$c Constraints for the relevant CYS-CYS
2006 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
2007 c$$$ f=f+tempf*tempf
2008 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
2010 c$$$c$$$ if (nf.ne.nfl) then
2011 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
2012 c$$$c$$$ + f,dist(5+nres,14+nres)
2020 c$$$C-----------------------------------------------------------------------------
2021 c$$$C-----------------------------------------------------------------------------
2022 subroutine triple_ssbond_ene(resi,resj,resk,eij)
2023 include 'DIMENSIONS'
2024 include 'COMMON.SBRIDGE'
2025 include 'COMMON.CHAIN'
2026 include 'COMMON.DERIV'
2027 include 'COMMON.LOCAL'
2028 include 'COMMON.INTERACT'
2029 include 'COMMON.VAR'
2030 include 'COMMON.IOUNITS'
2031 include 'COMMON.CALC'
2034 C include 'COMMON.MD'
2038 c External functions
2039 double precision h_base
2043 integer resi,resj,resk
2046 double precision eij,eij1,eij2,eij3
2050 c integer itypi,itypj,k,l
2051 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2052 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2053 double precision xik,yik,zik,xjk,yjk,zjk
2054 double precision sig0ij,ljd,sig,fac,e1,e2
2055 double precision dcosom1(3),dcosom2(3),ed
2056 double precision pom1,pom2
2057 double precision ljA,ljB,ljXs
2058 double precision d_ljB(1:3)
2059 double precision ssA,ssB,ssC,ssXs
2060 double precision ssxm,ljxm,ssm,ljm
2061 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2066 C write(iout,*) resi,resj,resk
2068 dxi=dc_norm(1,nres+i)
2069 dyi=dc_norm(2,nres+i)
2070 dzi=dc_norm(3,nres+i)
2071 dsci_inv=vbld_inv(i+nres)
2081 dxj=dc_norm(1,nres+j)
2082 dyj=dc_norm(2,nres+j)
2083 dzj=dc_norm(3,nres+j)
2084 dscj_inv=vbld_inv(j+nres)
2090 dxk=dc_norm(1,nres+k)
2091 dyk=dc_norm(2,nres+k)
2092 dzk=dc_norm(3,nres+k)
2093 dscj_inv=vbld_inv(k+nres)
2103 rrij=(xij*xij+yij*yij+zij*zij)
2104 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2105 rrik=(xik*xik+yik*yik+zik*zik)
2107 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2109 C there are three combination of distances for each trisulfide bonds
2110 C The first case the ith atom is the center
2111 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2112 C distance y is second distance the a,b,c,d are parameters derived for
2113 C this problem d parameter was set as a penalty currenlty set to 1.
2114 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2115 C second case jth atom is center
2116 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2117 C the third case kth atom is the center
2118 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2123 C write(iout,*)i,j,k,eij
2124 C The energy penalty calculated now time for the gradient part
2125 C derivative over rij
2126 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2127 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2132 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2133 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2136 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2137 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2139 C now derivative over rik
2140 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2141 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2146 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2147 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2150 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2151 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2153 C now derivative over rjk
2154 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2155 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2160 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2161 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2164 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2165 gvdwc(l,k)=gvdwc(l,k)+gg(l)