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)
155 xj=c(1,nres+j)-c(1,nres+i)
156 yj=c(2,nres+j)-c(2,nres+i)
157 zj=c(3,nres+j)-c(3,nres+i)
158 dxj=dc_norm(1,nres+j)
159 dyj=dc_norm(2,nres+j)
160 dzj=dc_norm(3,nres+j)
161 dscj_inv=vbld_inv(j+nres)
163 chi1=chi(itypi,itypj)
164 chi2=chi(itypj,itypi)
171 alf12=0.5D0*(alf1+alf2)
173 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
174 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
175 c The following are set in sc_angular
179 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
180 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
181 c om12=dxi*dxj+dyi*dyj+dzi*dzj
183 rij=1.0D0/rij ! Reset this so it makes sense
185 sig0ij=sigma(itypi,itypj)
186 sig=sig0ij*dsqrt(1.0D0/sigsq)
189 ljA=eps1*eps2rt**2*eps3rt**2
190 ljB=ljA*bb(itypi,itypj)
191 ljA=ljA*aa(itypi,itypj)
192 ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
197 deltat12=om2-om1+2.0d0
202 & +akth*(deltat1*deltat1+deltat2*deltat2)
203 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
204 ssxm=ssXs-0.5D0*ssB/ssA
207 c$$$c Some extra output
208 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
209 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
210 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
211 c$$$ if (ssx0.gt.0.0d0) then
212 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
216 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
217 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
218 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
220 c-------END TESTING CODE
223 c Stop and plot energy and derivative as a function of distance
225 ssm=ssC-0.25D0*ssB*ssB/ssA
226 ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
228 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
236 if (.not.checkstop) then
243 if (checkstop) rij=(ssxm-1.0d0)+
244 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
245 c-------END TESTING CODE
247 if (rij.gt.ljxm) then
250 fac=(1.0D0/ljd)**expon
251 e1=fac*fac*aa(itypi,itypj)
252 e2=fac*bb(itypi,itypj)
253 eij=eps1*eps2rt*eps3rt*(e1+e2)
254 C write(iout,*) eij,'TU?1'
257 eij=eij*eps2rt*eps3rt
260 e1=e1*eps1*eps2rt**2*eps3rt**2
261 ed=-expon*(e1+eij)/ljd
263 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
264 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
265 eom12=eij*eps1_om12+eps2der*eps2rt_om12
266 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
267 else if (rij.lt.ssxm) then
270 eij=ssA*ssd*ssd+ssB*ssd+ssC
271 C write(iout,*) 'TU?2',ssc,ssd
272 ed=2*akcm*ssd+akct*deltat12
274 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
275 eom1=-2*akth*deltat1-pom1-om2*pom2
276 eom2= 2*akth*deltat2+pom1-om1*pom2
279 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
281 d_ssxm(1)=0.5D0*akct/ssA
285 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
286 d_ljxm(2)=d_ljxm(1)*sigsq_om2
287 d_ljxm(3)=d_ljxm(1)*sigsq_om12
288 d_ljxm(1)=d_ljxm(1)*sigsq_om1
290 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
293 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
297 ssm=ssC-0.25D0*ssB*ssB/ssA
298 d_ssm(1)=0.5D0*akct*ssB/ssA
299 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
300 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
302 f1=(rij-xm)/(ssxm-xm)
303 f2=(rij-ssxm)/(xm-ssxm)
307 C write(iout,*) eij,'TU?3'
308 delta_inv=1.0d0/(xm-ssxm)
309 deltasq_inv=delta_inv*delta_inv
311 fac1=deltasq_inv*fac*(xm-rij)
312 fac2=deltasq_inv*fac*(rij-ssxm)
313 ed=delta_inv*(Ht*hd2-ssm*hd1)
314 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
315 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
316 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
319 ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
320 d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
321 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
322 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
324 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
325 f1=(rij-ljxm)/(xm-ljxm)
326 f2=(rij-xm)/(ljxm-xm)
330 C write(iout,*) 'TU?4',ssA
331 delta_inv=1.0d0/(ljxm-xm)
332 deltasq_inv=delta_inv*delta_inv
334 fac1=deltasq_inv*fac*(ljxm-rij)
335 fac2=deltasq_inv*fac*(rij-xm)
336 ed=delta_inv*(ljm*hd2-Ht*hd1)
337 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
338 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
339 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
341 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
343 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
349 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
350 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
351 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
353 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
354 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
355 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
356 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
359 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
361 c$$$ d_ljm(k)=ljm*d_ljB(k)
365 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
366 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
367 c$$$ d_ss(2)=akct*ssd
368 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
369 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
372 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
373 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
374 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
376 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
377 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
379 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
381 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
382 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
383 c$$$ h1=h_base(f1,hd1)
384 c$$$ h2=h_base(f2,hd2)
385 c$$$ eij=ss*h1+ljf*h2
386 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
387 c$$$ deltasq_inv=delta_inv*delta_inv
388 c$$$ fac=ljf*hd2-ss*hd1
389 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
390 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
391 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
392 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
393 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
394 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
395 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
397 c$$$ havebond=.false.
398 c$$$ if (ed.gt.0.0d0) havebond=.true.
399 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
402 write(iout,*) 'havebond',havebond
406 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
407 c write(iout,'(a15,f12.2,f8.1,2i5)')
408 c & "SSBOND_E_FORM",totT,t_bath,i,j
412 dyn_ssbond_ij(i,j)=eij
413 else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
414 dyn_ssbond_ij(i,j)=1.0d300
417 c write(iout,'(a15,f12.2,f8.1,2i5)')
418 c & "SSBOND_E_BREAK",totT,t_bath,i,j
425 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
426 & "CHECKSTOP",rij,eij,ed
431 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
438 c-------END TESTING CODE
441 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
442 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
445 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
448 gvdwx(k,i)=gvdwx(k,i)-gg(k)
449 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
450 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
451 gvdwx(k,j)=gvdwx(k,j)+gg(k)
452 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
453 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
457 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
462 gvdwc(l,i)=gvdwc(l,i)-gg(l)
463 gvdwc(l,j)=gvdwc(l,j)+gg(l)
469 C-----------------------------------------------------------------------------
471 double precision function h_base(x,deriv)
472 c A smooth function going 0->1 in range [0,1]
473 c It should NOT be called outside range [0,1], it will not work there.
480 double precision deriv
486 c Two parabolas put together. First derivative zero at extrema
487 c$$$ if (x.lt.0.5D0) then
488 c$$$ h_base=2.0D0*x*x
492 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
493 c$$$ deriv=4.0D0*deriv
496 c Third degree polynomial. First derivative zero at extrema
497 h_base=x*x*(3.0d0-2.0d0*x)
498 deriv=6.0d0*x*(1.0d0-x)
500 c Fifth degree polynomial. First and second derivatives zero at extrema
502 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
504 c$$$ deriv=deriv*deriv
505 c$$$ deriv=30.0d0*xsq*deriv
510 c----------------------------------------------------------------------------
512 subroutine dyn_set_nss
513 c Adjust nss and other relevant variables based on dyn_ssbond_ij
521 include 'COMMON.SBRIDGE'
522 include 'COMMON.CHAIN'
523 include 'COMMON.IOUNITS'
524 C include 'COMMON.SETUP'
527 C include 'COMMON.MD'
532 double precision emin
534 integer diff,allflag(maxdim),allnss,
535 & allihpb(maxdim),alljhpb(maxdim),
536 & newnss,newihpb(maxdim),newjhpb(maxdim)
538 integer i_newnss(1024),displ(0:1024)
539 integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
544 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
553 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
557 if (allflag(i).eq.0 .and.
558 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
559 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
563 if (emin.lt.1.0d300) then
566 if (allflag(i).eq.0 .and.
567 & (allihpb(i).eq.allihpb(imin) .or.
568 & alljhpb(i).eq.allihpb(imin) .or.
569 & allihpb(i).eq.alljhpb(imin) .or.
570 & alljhpb(i).eq.alljhpb(imin))) then
577 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
581 if (allflag(i).eq.1) then
583 newihpb(newnss)=allihpb(i)
584 newjhpb(newnss)=alljhpb(i)
589 if (nfgtasks.gt.1)then
591 call MPI_Reduce(newnss,g_newnss,1,
592 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
593 call MPI_Gather(newnss,1,MPI_INTEGER,
594 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
597 displ(i)=i_newnss(i-1)+displ(i-1)
599 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
600 & g_newihpb,i_newnss,displ,MPI_INTEGER,
602 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
603 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
605 if(fg_rank.eq.0) then
606 c print *,'g_newnss',g_newnss
607 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
608 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
611 newihpb(i)=g_newihpb(i)
612 newjhpb(i)=g_newjhpb(i)
620 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
625 if (idssb(i).eq.newihpb(j) .and.
626 & jdssb(i).eq.newjhpb(j)) found=.true.
630 c if (.not.found.and.fg_rank.eq.0)
631 c & write(iout,'(a15,f12.2,f8.1,2i5)')
632 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
640 if (newihpb(i).eq.idssb(j) .and.
641 & newjhpb(i).eq.jdssb(j)) found=.true.
645 c if (.not.found.and.fg_rank.eq.0)
646 c & write(iout,'(a15,f12.2,f8.1,2i5)')
647 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
661 c----------------------------------------------------------------------------
664 subroutine read_ssHist
669 include "DIMENSIONS.FREE"
670 include 'COMMON.FREE'
674 character*80 controlcard
677 call card_concat(controlcard,.true.)
679 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
686 c----------------------------------------------------------------------------
689 C-----------------------------------------------------------------------------
690 C-----------------------------------------------------------------------------
691 C-----------------------------------------------------------------------------
692 C-----------------------------------------------------------------------------
693 C-----------------------------------------------------------------------------
694 C-----------------------------------------------------------------------------
695 C-----------------------------------------------------------------------------
697 c$$$c-----------------------------------------------------------------------------
699 c$$$ subroutine ss_relax(i_in,j_in)
703 c$$$ include 'DIMENSIONS'
704 c$$$ include 'COMMON.VAR'
705 c$$$ include 'COMMON.CHAIN'
706 c$$$ include 'COMMON.IOUNITS'
707 c$$$ include 'COMMON.INTERACT'
709 c$$$c Input arguments
710 c$$$ integer i_in,j_in
712 c$$$c Local variables
713 c$$$ integer i,iretcode,nfun_sc
715 c$$$ double precision var(maxvar),e_sc,etot
722 c$$$ mask_side(i_in)=1
723 c$$$ mask_side(j_in)=1
725 c$$$c Minimize the two selected side-chains
726 c$$$ call overlap_sc(scfail) ! Better not fail!
727 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
734 c$$$c-------------------------------------------------------------
736 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
737 c$$$c Minimize side-chains only, starting from geom but without modifying
739 c$$$c If mask_r is already set, only the selected side-chains are minimized,
740 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
744 c$$$ include 'DIMENSIONS'
745 c$$$ include 'COMMON.IOUNITS'
746 c$$$ include 'COMMON.VAR'
747 c$$$ include 'COMMON.CHAIN'
748 c$$$ include 'COMMON.GEO'
749 c$$$ include 'COMMON.MINIM'
751 c$$$ common /srutu/ icall
753 c$$$c Output arguments
754 c$$$ double precision etot_sc
755 c$$$ integer iretcode,nfun
757 c$$$c External functions/subroutines
758 c$$$ external func_sc,grad_sc,fdum
760 c$$$c Local variables
762 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
764 c$$$ double precision rdum(1)
765 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
767 c$$$ integer i,nvar_restr
770 c$$$cmc start_minim=.true.
771 c$$$ call deflt(2,iv,liv,lv,v)
772 c$$$* 12 means fresh start, dont call deflt
774 c$$$* max num of fun calls
775 c$$$ if (maxfun.eq.0) maxfun=500
777 c$$$* max num of iterations
778 c$$$ if (maxmin.eq.0) maxmin=1000
780 c$$$* controls output
782 c$$$* selects output unit
784 c$$$c iv(21)=iout ! DEBUG
785 c$$$c iv(21)=8 ! DEBUG
786 c$$$* 1 means to print out result
788 c$$$c iv(22)=1 ! DEBUG
789 c$$$* 1 means to print out summary stats
791 c$$$c iv(23)=1 ! DEBUG
792 c$$$* 1 means to print initial x and d
794 c$$$c iv(24)=1 ! DEBUG
795 c$$$* min val for v(radfac) default is 0.1
797 c$$$* max val for v(radfac) default is 4.0
800 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
801 c$$$* the sumsl default is 0.1
803 c$$$* false conv if (act fnctn decrease) .lt. v(34)
804 c$$$* the sumsl default is 100*machep
805 c$$$ v(34)=v(34)/100.0D0
806 c$$$* absolute convergence
807 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
809 c$$$* relative convergence
810 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
812 c$$$* controls initial step size
814 c$$$* large vals of d correspond to small components of step
818 c$$$ do i=nphi+1,nvar
822 c$$$ call geom_to_var(nvar,x)
823 c$$$ IF (mask_r) THEN
824 c$$$ do i=1,nres ! Just in case...
828 c$$$ call x2xx(x,xx,nvar_restr)
829 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
830 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
833 c$$$c When minimizing ALL side-chains, etotal_sc is a little
834 c$$$c faster if we don't set mask_r
840 c$$$ call x2xx(x,xx,nvar_restr)
841 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
842 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
845 c$$$ call var_to_geom(nvar,x)
846 c$$$ call chainbuild_sc
853 c$$$C--------------------------------------------------------------------------
855 c$$$ subroutine chainbuild_sc
857 c$$$ include 'DIMENSIONS'
858 c$$$ include 'COMMON.VAR'
859 c$$$ include 'COMMON.INTERACT'
861 c$$$c Local variables
866 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
867 c$$$ call locate_side_chain(i)
874 c$$$C--------------------------------------------------------------------------
876 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
880 c$$$ include 'DIMENSIONS'
881 c$$$ include 'COMMON.DERIV'
882 c$$$ include 'COMMON.VAR'
883 c$$$ include 'COMMON.MINIM'
884 c$$$ include 'COMMON.IOUNITS'
886 c$$$c Input arguments
888 c$$$ double precision x(maxvar)
889 c$$$ double precision ufparm
892 c$$$c Input/Output arguments
894 c$$$ integer uiparm(1)
895 c$$$ double precision urparm(1)
897 c$$$c Output arguments
898 c$$$ double precision f
900 c$$$c Local variables
901 c$$$ double precision energia(0:n_ene)
903 c$$$c Variables used to intercept NaNs
904 c$$$ double precision x_sum
913 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
916 c$$$ x_sum=x_sum+x(i_NAN)
918 c$$$c Calculate the energy only if the coordinates are ok
919 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
920 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
926 c$$$ call var_to_geom_restr(n,x)
928 c$$$ call chainbuild_sc
929 c$$$ call etotal_sc(energia(0))
931 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
940 c$$$c-------------------------------------------------------
942 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
946 c$$$ include 'DIMENSIONS'
947 c$$$ include 'COMMON.CHAIN'
948 c$$$ include 'COMMON.DERIV'
949 c$$$ include 'COMMON.VAR'
950 c$$$ include 'COMMON.INTERACT'
951 c$$$ include 'COMMON.MINIM'
953 c$$$c Input arguments
955 c$$$ double precision x(maxvar)
956 c$$$ double precision ufparm
959 c$$$c Input/Output arguments
961 c$$$ integer uiparm(1)
962 c$$$ double precision urparm(1)
964 c$$$c Output arguments
965 c$$$ double precision g(maxvar)
967 c$$$c Local variables
968 c$$$ double precision f,gphii,gthetai,galphai,gomegai
969 c$$$ integer ig,ind,i,j,k,igall,ij
973 c$$$ if (nf-nfl+1) 20,30,40
974 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
975 c$$$c write (iout,*) 'grad 20'
976 c$$$ if (nf.eq.0) return
978 c$$$ 30 call var_to_geom_restr(n,x)
979 c$$$ call chainbuild_sc
981 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
985 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
991 c$$$ IF (mask_phi(i+2).eq.1) THEN
996 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
997 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1003 c$$$ ind=ind+nres-1-i
1010 c$$$ IF (mask_theta(i+2).eq.1) THEN
1013 c$$$ do j=i+1,nres-1
1016 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1017 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1022 c$$$ ind=ind+nres-1-i
1027 c$$$ if (itype(i).ne.10) then
1028 c$$$ IF (mask_side(i).eq.1) THEN
1032 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1041 c$$$ if (itype(i).ne.10) then
1042 c$$$ IF (mask_side(i).eq.1) THEN
1046 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1054 c$$$C Add the components corresponding to local energy terms.
1061 c$$$ if (mask_phi(i).eq.1) then
1063 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1069 c$$$ if (mask_theta(i).eq.1) then
1071 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1077 c$$$ if (itype(i).ne.10) then
1079 c$$$ if (mask_side(i).eq.1) then
1081 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1088 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1094 c$$$C-----------------------------------------------------------------------------
1096 c$$$ subroutine etotal_sc(energy_sc)
1100 c$$$ include 'DIMENSIONS'
1101 c$$$ include 'COMMON.VAR'
1102 c$$$ include 'COMMON.INTERACT'
1103 c$$$ include 'COMMON.DERIV'
1104 c$$$ include 'COMMON.FFIELD'
1106 c$$$c Output arguments
1107 c$$$ double precision energy_sc(0:n_ene)
1109 c$$$c Local variables
1110 c$$$ double precision evdw,escloc
1115 c$$$ energy_sc(i)=0.0D0
1118 c$$$ if (mask_r) then
1119 c$$$ call egb_sc(evdw)
1120 c$$$ call esc_sc(escloc)
1123 c$$$ call esc(escloc)
1126 c$$$ if (evdw.eq.1.0D20) then
1127 c$$$ energy_sc(0)=evdw
1129 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1131 c$$$ energy_sc(1)=evdw
1132 c$$$ energy_sc(12)=escloc
1135 c$$$C Sum up the components of the Cartesian gradient.
1139 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1146 c$$$C-----------------------------------------------------------------------------
1148 c$$$ subroutine egb_sc(evdw)
1150 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1151 c$$$C assuming the Gay-Berne potential of interaction.
1153 c$$$ implicit real*8 (a-h,o-z)
1154 c$$$ include 'DIMENSIONS'
1155 c$$$ include 'COMMON.GEO'
1156 c$$$ include 'COMMON.VAR'
1157 c$$$ include 'COMMON.LOCAL'
1158 c$$$ include 'COMMON.CHAIN'
1159 c$$$ include 'COMMON.DERIV'
1160 c$$$ include 'COMMON.NAMES'
1161 c$$$ include 'COMMON.INTERACT'
1162 c$$$ include 'COMMON.IOUNITS'
1163 c$$$ include 'COMMON.CALC'
1164 c$$$ include 'COMMON.CONTROL'
1167 c$$$ energy_dec=.false.
1168 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1171 c$$$c if (icall.eq.0) lprn=.false.
1173 c$$$ do i=iatsc_s,iatsc_e
1175 c$$$ itypi1=itype(i+1)
1179 c$$$ dxi=dc_norm(1,nres+i)
1180 c$$$ dyi=dc_norm(2,nres+i)
1181 c$$$ dzi=dc_norm(3,nres+i)
1182 c$$$c dsci_inv=dsc_inv(itypi)
1183 c$$$ dsci_inv=vbld_inv(i+nres)
1184 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1185 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1187 c$$$C Calculate SC interaction energy.
1189 c$$$ do iint=1,nint_gr(i)
1190 c$$$ do j=istart(i,iint),iend(i,iint)
1191 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1194 c$$$c dscj_inv=dsc_inv(itypj)
1195 c$$$ dscj_inv=vbld_inv(j+nres)
1196 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1197 c$$$c & 1.0d0/vbld(j+nres)
1198 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1199 c$$$ sig0ij=sigma(itypi,itypj)
1200 c$$$ chi1=chi(itypi,itypj)
1201 c$$$ chi2=chi(itypj,itypi)
1202 c$$$ chi12=chi1*chi2
1203 c$$$ chip1=chip(itypi)
1204 c$$$ chip2=chip(itypj)
1205 c$$$ chip12=chip1*chip2
1206 c$$$ alf1=alp(itypi)
1207 c$$$ alf2=alp(itypj)
1208 c$$$ alf12=0.5D0*(alf1+alf2)
1209 c$$$C For diagnostics only!!!
1219 c$$$ xj=c(1,nres+j)-xi
1220 c$$$ yj=c(2,nres+j)-yi
1221 c$$$ zj=c(3,nres+j)-zi
1222 c$$$ dxj=dc_norm(1,nres+j)
1223 c$$$ dyj=dc_norm(2,nres+j)
1224 c$$$ dzj=dc_norm(3,nres+j)
1225 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1226 c$$$c write (iout,*) "j",j," dc_norm",
1227 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1228 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1229 c$$$ rij=dsqrt(rrij)
1230 c$$$C Calculate angle-dependent terms of energy and contributions to their
1232 c$$$ call sc_angular
1233 c$$$ sigsq=1.0D0/sigsq
1234 c$$$ sig=sig0ij*dsqrt(sigsq)
1235 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1236 c$$$c for diagnostics; uncomment
1237 c$$$c rij_shift=1.2*sig0ij
1238 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1239 c$$$ if (rij_shift.le.0.0D0) then
1241 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1242 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1243 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1246 c$$$ sigder=-sig*sigsq
1247 c$$$c---------------------------------------------------------------
1248 c$$$ rij_shift=1.0D0/rij_shift
1249 c$$$ fac=rij_shift**expon
1250 c$$$ e1=fac*fac*aa(itypi,itypj)
1251 c$$$ e2=fac*bb(itypi,itypj)
1252 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1253 c$$$ eps2der=evdwij*eps3rt
1254 c$$$ eps3der=evdwij*eps2rt
1255 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1256 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1257 c$$$ evdwij=evdwij*eps2rt*eps3rt
1258 c$$$ evdw=evdw+evdwij
1260 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1261 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1262 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1263 c$$$ & restyp(itypi),i,restyp(itypj),j,
1264 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1265 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1266 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1270 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1271 c$$$ & 'evdw',i,j,evdwij
1273 c$$$C Calculate gradient components.
1274 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1275 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1276 c$$$ sigder=fac*sigder
1279 c$$$C Calculate the radial part of the gradient
1283 c$$$C Calculate angular part of the gradient.
1289 c$$$ energy_dec=.false.
1293 c$$$c-----------------------------------------------------------------------------
1295 c$$$ subroutine esc_sc(escloc)
1296 c$$$C Calculate the local energy of a side chain and its derivatives in the
1297 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1298 c$$$C ALPHA and OMEGA.
1299 c$$$ implicit real*8 (a-h,o-z)
1300 c$$$ include 'DIMENSIONS'
1301 c$$$ include 'COMMON.GEO'
1302 c$$$ include 'COMMON.LOCAL'
1303 c$$$ include 'COMMON.VAR'
1304 c$$$ include 'COMMON.INTERACT'
1305 c$$$ include 'COMMON.DERIV'
1306 c$$$ include 'COMMON.CHAIN'
1307 c$$$ include 'COMMON.IOUNITS'
1308 c$$$ include 'COMMON.NAMES'
1309 c$$$ include 'COMMON.FFIELD'
1310 c$$$ include 'COMMON.CONTROL'
1311 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1312 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1313 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1314 c$$$ delta=0.02d0*pi
1316 c$$$c write (iout,'(a)') 'ESC'
1317 c$$$ do i=loc_start,loc_end
1318 c$$$ IF (mask_side(i).eq.1) THEN
1320 c$$$ if (it.eq.10) goto 1
1321 c$$$ nlobit=nlob(it)
1322 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1323 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1324 c$$$ theti=theta(i+1)-pipol
1325 c$$$ x(1)=dtan(theti)
1329 c$$$ if (x(2).gt.pi-delta) then
1331 c$$$ xtemp(2)=pi-delta
1333 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1335 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1336 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1337 c$$$ & escloci,dersc(2))
1338 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1339 c$$$ & ddersc0(1),dersc(1))
1340 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1341 c$$$ & ddersc0(3),dersc(3))
1342 c$$$ xtemp(2)=pi-delta
1343 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1345 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1346 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1347 c$$$ & dersc0(2),esclocbi,dersc02)
1348 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1349 c$$$ & dersc12,dersc01)
1350 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1351 c$$$ dersc0(1)=dersc01
1352 c$$$ dersc0(2)=dersc02
1353 c$$$ dersc0(3)=0.0d0
1355 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1357 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1358 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1359 c$$$c & esclocbi,ss,ssd
1360 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1361 c$$$c escloci=esclocbi
1362 c$$$c write (iout,*) escloci
1363 c$$$ else if (x(2).lt.delta) then
1367 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1369 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1370 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1371 c$$$ & escloci,dersc(2))
1372 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1373 c$$$ & ddersc0(1),dersc(1))
1374 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1375 c$$$ & ddersc0(3),dersc(3))
1377 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1379 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1380 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1381 c$$$ & dersc0(2),esclocbi,dersc02)
1382 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1383 c$$$ & dersc12,dersc01)
1384 c$$$ dersc0(1)=dersc01
1385 c$$$ dersc0(2)=dersc02
1386 c$$$ dersc0(3)=0.0d0
1387 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1389 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1391 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1392 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1393 c$$$c & esclocbi,ss,ssd
1394 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1395 c$$$c write (iout,*) escloci
1397 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1400 c$$$ escloc=escloc+escloci
1401 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1402 c$$$ & 'escloc',i,escloci
1403 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1405 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1406 c$$$ & wscloc*dersc(1)
1407 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1408 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1415 c$$$C-----------------------------------------------------------------------------
1417 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1419 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1420 c$$$C assuming the Gay-Berne potential of interaction.
1422 c$$$ implicit real*8 (a-h,o-z)
1423 c$$$ include 'DIMENSIONS'
1424 c$$$ include 'COMMON.GEO'
1425 c$$$ include 'COMMON.VAR'
1426 c$$$ include 'COMMON.LOCAL'
1427 c$$$ include 'COMMON.CHAIN'
1428 c$$$ include 'COMMON.DERIV'
1429 c$$$ include 'COMMON.NAMES'
1430 c$$$ include 'COMMON.INTERACT'
1431 c$$$ include 'COMMON.IOUNITS'
1432 c$$$ include 'COMMON.CALC'
1433 c$$$ include 'COMMON.CONTROL'
1436 c$$$ energy_dec=.false.
1437 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1441 c$$$c$$$ do i=iatsc_s,iatsc_e
1444 c$$$ itypi1=itype(i+1)
1448 c$$$ dxi=dc_norm(1,nres+i)
1449 c$$$ dyi=dc_norm(2,nres+i)
1450 c$$$ dzi=dc_norm(3,nres+i)
1451 c$$$c dsci_inv=dsc_inv(itypi)
1452 c$$$ dsci_inv=vbld_inv(i+nres)
1453 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1454 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1456 c$$$C Calculate SC interaction energy.
1458 c$$$c$$$ do iint=1,nint_gr(i)
1459 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1463 c$$$c dscj_inv=dsc_inv(itypj)
1464 c$$$ dscj_inv=vbld_inv(j+nres)
1465 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1466 c$$$c & 1.0d0/vbld(j+nres)
1467 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1468 c$$$ sig0ij=sigma(itypi,itypj)
1469 c$$$ chi1=chi(itypi,itypj)
1470 c$$$ chi2=chi(itypj,itypi)
1471 c$$$ chi12=chi1*chi2
1472 c$$$ chip1=chip(itypi)
1473 c$$$ chip2=chip(itypj)
1474 c$$$ chip12=chip1*chip2
1475 c$$$ alf1=alp(itypi)
1476 c$$$ alf2=alp(itypj)
1477 c$$$ alf12=0.5D0*(alf1+alf2)
1478 c$$$C For diagnostics only!!!
1488 c$$$ xj=c(1,nres+j)-xi
1489 c$$$ yj=c(2,nres+j)-yi
1490 c$$$ zj=c(3,nres+j)-zi
1491 c$$$ dxj=dc_norm(1,nres+j)
1492 c$$$ dyj=dc_norm(2,nres+j)
1493 c$$$ dzj=dc_norm(3,nres+j)
1494 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1495 c$$$c write (iout,*) "j",j," dc_norm",
1496 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1497 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1498 c$$$ rij=dsqrt(rrij)
1499 c$$$C Calculate angle-dependent terms of energy and contributions to their
1501 c$$$ call sc_angular
1502 c$$$ sigsq=1.0D0/sigsq
1503 c$$$ sig=sig0ij*dsqrt(sigsq)
1504 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1505 c$$$c for diagnostics; uncomment
1506 c$$$c rij_shift=1.2*sig0ij
1507 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1508 c$$$ if (rij_shift.le.0.0D0) then
1510 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1511 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1512 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1515 c$$$ sigder=-sig*sigsq
1516 c$$$c---------------------------------------------------------------
1517 c$$$ rij_shift=1.0D0/rij_shift
1518 c$$$ fac=rij_shift**expon
1519 c$$$ e1=fac*fac*aa(itypi,itypj)
1520 c$$$ e2=fac*bb(itypi,itypj)
1521 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1522 c$$$ eps2der=evdwij*eps3rt
1523 c$$$ eps3der=evdwij*eps2rt
1524 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1525 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1526 c$$$ evdwij=evdwij*eps2rt*eps3rt
1527 c$$$ evdw=evdw+evdwij
1529 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1530 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1531 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1532 c$$$ & restyp(itypi),i,restyp(itypj),j,
1533 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1534 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1535 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1539 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1540 c$$$ & 'evdw',i,j,evdwij
1542 c$$$C Calculate gradient components.
1543 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1544 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1545 c$$$ sigder=fac*sigder
1548 c$$$C Calculate the radial part of the gradient
1552 c$$$C Calculate angular part of the gradient.
1555 c$$$c$$$ enddo ! iint
1557 c$$$ energy_dec=.false.
1561 c$$$C-----------------------------------------------------------------------------
1563 c$$$ subroutine perturb_side_chain(i,angle)
1567 c$$$ include 'DIMENSIONS'
1568 c$$$ include 'COMMON.CHAIN'
1569 c$$$ include 'COMMON.GEO'
1570 c$$$ include 'COMMON.VAR'
1571 c$$$ include 'COMMON.LOCAL'
1572 c$$$ include 'COMMON.IOUNITS'
1574 c$$$c External functions
1575 c$$$ external ran_number
1576 c$$$ double precision ran_number
1578 c$$$c Input arguments
1580 c$$$ double precision angle ! In degrees
1582 c$$$c Local variables
1584 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1588 c$$$ rad_ang=angle*deg2rad
1591 c$$$ do while (length.lt.0.01)
1592 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1593 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1594 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1595 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1596 c$$$ + rand_v(3)*rand_v(3)
1597 c$$$ length=sqrt(length)
1598 c$$$ rand_v(1)=rand_v(1)/length
1599 c$$$ rand_v(2)=rand_v(2)/length
1600 c$$$ rand_v(3)=rand_v(3)/length
1601 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1602 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1603 c$$$ length=1.0D0-cost*cost
1604 c$$$ if (length.lt.0.0D0) length=0.0D0
1605 c$$$ length=sqrt(length)
1606 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1607 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1608 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1610 c$$$ rand_v(1)=rand_v(1)/length
1611 c$$$ rand_v(2)=rand_v(2)/length
1612 c$$$ rand_v(3)=rand_v(3)/length
1614 c$$$ cost=dcos(rad_ang)
1615 c$$$ sint=dsin(rad_ang)
1616 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1617 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1618 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1619 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1620 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1621 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1622 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1623 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1624 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1626 c$$$ call chainbuild_cart
1631 c$$$c----------------------------------------------------------------------------
1633 c$$$ subroutine ss_relax3(i_in,j_in)
1637 c$$$ include 'DIMENSIONS'
1638 c$$$ include 'COMMON.VAR'
1639 c$$$ include 'COMMON.CHAIN'
1640 c$$$ include 'COMMON.IOUNITS'
1641 c$$$ include 'COMMON.INTERACT'
1643 c$$$c External functions
1644 c$$$ external ran_number
1645 c$$$ double precision ran_number
1647 c$$$c Input arguments
1648 c$$$ integer i_in,j_in
1650 c$$$c Local variables
1651 c$$$ double precision energy_sc(0:n_ene),etot
1652 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1653 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1654 c$$$ integer n,i_pert,i
1655 c$$$ logical notdone
1664 c$$$ mask_side(i_in)=1
1665 c$$$ mask_side(j_in)=1
1667 c$$$ call etotal_sc(energy_sc)
1668 c$$$ etot=energy_sc(0)
1669 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1670 c$$$c + energy_sc(1),energy_sc(12)
1674 c$$$ do while (notdone)
1675 c$$$ if (mod(n,2).eq.0) then
1683 c$$$ org_dc(i)=dc(i,i_pert+nres)
1684 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1685 c$$$ org_c(i)=c(i,i_pert+nres)
1687 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1688 c$$$ call perturb_side_chain(i_pert,ang_pert)
1689 c$$$ call etotal_sc(energy_sc)
1690 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1691 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1692 c$$$ if (rand_fact.lt.exp_fact) then
1693 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1694 c$$$c + energy_sc(1),energy_sc(12)
1695 c$$$ etot=energy_sc(0)
1697 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1698 c$$$c + energy_sc(1),energy_sc(12)
1700 c$$$ dc(i,i_pert+nres)=org_dc(i)
1701 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1702 c$$$ c(i,i_pert+nres)=org_c(i)
1706 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1714 c$$$c----------------------------------------------------------------------------
1716 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1718 c$$$ include 'DIMENSIONS'
1720 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1721 c$$$*********************************************************************
1722 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1723 c$$$* the calling subprogram. *
1724 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1725 c$$$* calculated in the usual pythagorean way. *
1726 c$$$* absolute convergence occurs when the function is within v(31) of *
1727 c$$$* zero. unless you know the minimum value in advance, abs convg *
1728 c$$$* is probably not useful. *
1729 c$$$* relative convergence is when the model predicts that the function *
1730 c$$$* will decrease by less than v(32)*abs(fun). *
1731 c$$$*********************************************************************
1732 c$$$ include 'COMMON.IOUNITS'
1733 c$$$ include 'COMMON.VAR'
1734 c$$$ include 'COMMON.GEO'
1735 c$$$ include 'COMMON.MINIM'
1736 c$$$ include 'COMMON.CHAIN'
1738 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1739 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1740 c$$$ + orig_ss_dist(maxres2,maxres2)
1742 c$$$ double precision etot
1743 c$$$ integer iretcode,nfun,i_in,j_in
1746 c$$$ double precision dist
1747 c$$$ external ss_func,fdum
1748 c$$$ double precision ss_func,fdum
1750 c$$$ integer iv(liv),uiparm(2)
1751 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1755 c$$$ call deflt(2,iv,liv,lv,v)
1756 c$$$* 12 means fresh start, dont call deflt
1758 c$$$* max num of fun calls
1759 c$$$ if (maxfun.eq.0) maxfun=500
1761 c$$$* max num of iterations
1762 c$$$ if (maxmin.eq.0) maxmin=1000
1764 c$$$* controls output
1766 c$$$* selects output unit
1769 c$$$* 1 means to print out result
1771 c$$$* 1 means to print out summary stats
1773 c$$$* 1 means to print initial x and d
1775 c$$$* min val for v(radfac) default is 0.1
1777 c$$$* max val for v(radfac) default is 4.0
1780 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1781 c$$$* the sumsl default is 0.1
1783 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1784 c$$$* the sumsl default is 100*machep
1785 c$$$ v(34)=v(34)/100.0D0
1786 c$$$* absolute convergence
1787 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1790 c$$$* relative convergence
1791 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1794 c$$$* controls initial step size
1796 c$$$* large vals of d correspond to small components of step
1803 c$$$ orig_ss_dc(j,i)=dc(j,i)
1806 c$$$ call geom_to_var(nvar,orig_ss_var)
1810 c$$$ orig_ss_dist(j,i)=dist(j,i)
1811 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1812 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1813 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1825 c$$$ if (ialph(i,1).gt.0) then
1828 c$$$ x(k)=dc(j,i+nres)
1835 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1838 c$$$ nfun=iv(6)+iv(30)
1848 c$$$ if (ialph(i,1).gt.0) then
1851 c$$$ dc(j,i+nres)=x(k)
1855 c$$$ call chainbuild_cart
1860 c$$$C-----------------------------------------------------------------------------
1862 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1864 c$$$ include 'DIMENSIONS'
1865 c$$$ include 'COMMON.DERIV'
1866 c$$$ include 'COMMON.IOUNITS'
1867 c$$$ include 'COMMON.VAR'
1868 c$$$ include 'COMMON.CHAIN'
1869 c$$$ include 'COMMON.INTERACT'
1870 c$$$ include 'COMMON.SBRIDGE'
1872 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1873 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1874 c$$$ + orig_ss_dist(maxres2,maxres2)
1877 c$$$ double precision x(maxres6)
1879 c$$$ double precision f
1880 c$$$ integer uiparm(2)
1881 c$$$ real*8 urparm(1)
1882 c$$$ external ufparm
1883 c$$$ double precision ufparm
1886 c$$$ double precision dist
1888 c$$$ integer i,j,k,ss_i,ss_j
1889 c$$$ double precision tempf,var(maxvar)
1904 c$$$ if (ialph(i,1).gt.0) then
1907 c$$$ dc(j,i+nres)=x(k)
1911 c$$$ call chainbuild_cart
1913 c$$$ call geom_to_var(nvar,var)
1915 c$$$c Constraints on all angles
1917 c$$$ tempf=var(i)-orig_ss_var(i)
1918 c$$$ f=f+tempf*tempf
1921 c$$$c Constraints on all distances
1923 c$$$ if (i.gt.1) then
1924 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1925 c$$$ f=f+tempf*tempf
1928 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1929 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1930 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1931 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1932 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1933 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1934 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1935 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1939 c$$$c Constraints for the relevant CYS-CYS
1940 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1941 c$$$ f=f+tempf*tempf
1942 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
1944 c$$$c$$$ if (nf.ne.nfl) then
1945 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1946 c$$$c$$$ + f,dist(5+nres,14+nres)
1954 c$$$C-----------------------------------------------------------------------------
1955 c$$$C-----------------------------------------------------------------------------
1956 subroutine triple_ssbond_ene(resi,resj,resk,eij)
1957 include 'DIMENSIONS'
1958 include 'COMMON.SBRIDGE'
1959 include 'COMMON.CHAIN'
1960 include 'COMMON.DERIV'
1961 include 'COMMON.LOCAL'
1962 include 'COMMON.INTERACT'
1963 include 'COMMON.VAR'
1964 include 'COMMON.IOUNITS'
1965 include 'COMMON.CALC'
1968 C include 'COMMON.MD'
1972 c External functions
1973 double precision h_base
1977 integer resi,resj,resk
1980 double precision eij,eij1,eij2,eij3
1984 c integer itypi,itypj,k,l
1985 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
1986 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
1987 double precision xik,yik,zik,xjk,yjk,zjk
1988 double precision sig0ij,ljd,sig,fac,e1,e2
1989 double precision dcosom1(3),dcosom2(3),ed
1990 double precision pom1,pom2
1991 double precision ljA,ljB,ljXs
1992 double precision d_ljB(1:3)
1993 double precision ssA,ssB,ssC,ssXs
1994 double precision ssxm,ljxm,ssm,ljm
1995 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2000 C write(iout,*) resi,resj,resk
2002 dxi=dc_norm(1,nres+i)
2003 dyi=dc_norm(2,nres+i)
2004 dzi=dc_norm(3,nres+i)
2005 dsci_inv=vbld_inv(i+nres)
2015 dxj=dc_norm(1,nres+j)
2016 dyj=dc_norm(2,nres+j)
2017 dzj=dc_norm(3,nres+j)
2018 dscj_inv=vbld_inv(j+nres)
2024 dxk=dc_norm(1,nres+k)
2025 dyk=dc_norm(2,nres+k)
2026 dzk=dc_norm(3,nres+k)
2027 dscj_inv=vbld_inv(k+nres)
2037 rrij=(xij*xij+yij*yij+zij*zij)
2038 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2039 rrik=(xik*xik+yik*yik+zik*zik)
2041 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2043 C there are three combination of distances for each trisulfide bonds
2044 C The first case the ith atom is the center
2045 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2046 C distance y is second distance the a,b,c,d are parameters derived for
2047 C this problem d parameter was set as a penalty currenlty set to 1.
2048 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2049 C second case jth atom is center
2050 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2051 C the third case kth atom is the center
2052 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2057 C write(iout,*)i,j,k,eij
2058 C The energy penalty calculated now time for the gradient part
2059 C derivative over rij
2060 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2061 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2066 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2067 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2070 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2071 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2073 C now derivative over rik
2074 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2075 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2080 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2081 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2084 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2085 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2087 C now derivative over rjk
2088 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2089 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2094 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2095 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2098 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2099 gvdwc(l,k)=gvdwc(l,k)+gg(l)