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'
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_aq(itypi,itypj)
191 ljA=ljA*aa_aq(itypi,itypj)
192 ljxm=ljXs+(-2.0D0*aa_aq(itypi,itypj)/
193 & bb_aq(itypi,itypj))**(1.0D0/6.0D0)
198 deltat12=om2-om1+2.0d0
203 & +akth*(deltat1*deltat1+deltat2*deltat2)
204 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
205 ssxm=ssXs-0.5D0*ssB/ssA
208 c$$$c Some extra output
209 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
210 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
211 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
212 c$$$ if (ssx0.gt.0.0d0) then
213 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
217 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
218 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
219 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
221 c-------END TESTING CODE
224 c Stop and plot energy and derivative as a function of distance
226 ssm=ssC-0.25D0*ssB*ssB/ssA
227 ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)
229 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
237 if (.not.checkstop) then
244 if (checkstop) rij=(ssxm-1.0d0)+
245 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
246 c-------END TESTING CODE
248 if (rij.gt.ljxm) then
251 fac=(1.0D0/ljd)**expon
252 e1=fac*fac*aa_aq(itypi,itypj)
253 e2=fac*bb_aq(itypi,itypj)
254 eij=eps1*eps2rt*eps3rt*(e1+e2)
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
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 delta_inv=1.0d0/(xm-ssxm)
308 deltasq_inv=delta_inv*delta_inv
310 fac1=deltasq_inv*fac*(xm-rij)
311 fac2=deltasq_inv*fac*(rij-ssxm)
312 ed=delta_inv*(Ht*hd2-ssm*hd1)
313 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
314 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
315 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
318 ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)
319 d_ljm(1)=-0.5D0*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)*ljB
320 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
321 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
323 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
324 f1=(rij-ljxm)/(xm-ljxm)
325 f2=(rij-xm)/(ljxm-xm)
329 delta_inv=1.0d0/(ljxm-xm)
330 deltasq_inv=delta_inv*delta_inv
332 fac1=deltasq_inv*fac*(ljxm-rij)
333 fac2=deltasq_inv*fac*(rij-xm)
334 ed=delta_inv*(ljm*hd2-Ht*hd1)
335 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
336 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
337 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
339 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
341 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
347 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
348 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
349 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
351 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
352 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
353 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
354 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
357 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
359 c$$$ d_ljm(k)=ljm*d_ljB(k)
363 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
364 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
365 c$$$ d_ss(2)=akct*ssd
366 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
367 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
370 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
371 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
372 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
374 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
375 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
377 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
379 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
380 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
381 c$$$ h1=h_base(f1,hd1)
382 c$$$ h2=h_base(f2,hd2)
383 c$$$ eij=ss*h1+ljf*h2
384 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
385 c$$$ deltasq_inv=delta_inv*delta_inv
386 c$$$ fac=ljf*hd2-ss*hd1
387 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
388 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
389 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
390 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
391 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
392 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
393 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
395 c$$$ havebond=.false.
396 c$$$ if (ed.gt.0.0d0) havebond=.true.
397 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
404 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
405 c write(iout,'(a15,f12.2,f8.1,2i5)')
406 c & "SSBOND_E_FORM",totT,t_bath,i,j
410 dyn_ssbond_ij(i,j)=eij
411 else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
412 dyn_ssbond_ij(i,j)=1.0d300
415 c write(iout,'(a15,f12.2,f8.1,2i5)')
416 c & "SSBOND_E_BREAK",totT,t_bath,i,j
423 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
424 & "CHECKSTOP",rij,eij,ed
429 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
436 c-------END TESTING CODE
439 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
440 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
443 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
446 gvdwx(k,i)=gvdwx(k,i)-gg(k)
447 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
448 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
449 gvdwx(k,j)=gvdwx(k,j)+gg(k)
450 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
451 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
455 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
460 gvdwc(l,i)=gvdwc(l,i)-gg(l)
461 gvdwc(l,j)=gvdwc(l,j)+gg(l)
467 C-----------------------------------------------------------------------------
469 double precision function h_base(x,deriv)
470 c A smooth function going 0->1 in range [0,1]
471 c It should NOT be called outside range [0,1], it will not work there.
478 double precision deriv
484 c Two parabolas put together. First derivative zero at extrema
485 c$$$ if (x.lt.0.5D0) then
486 c$$$ h_base=2.0D0*x*x
490 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
491 c$$$ deriv=4.0D0*deriv
494 c Third degree polynomial. First derivative zero at extrema
495 h_base=x*x*(3.0d0-2.0d0*x)
496 deriv=6.0d0*x*(1.0d0-x)
498 c Fifth degree polynomial. First and second derivatives zero at extrema
500 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
502 c$$$ deriv=deriv*deriv
503 c$$$ deriv=30.0d0*xsq*deriv
508 c----------------------------------------------------------------------------
510 subroutine dyn_set_nss
511 c Adjust nss and other relevant variables based on dyn_ssbond_ij
519 include 'COMMON.SBRIDGE'
520 include 'COMMON.CHAIN'
521 include 'COMMON.IOUNITS'
522 include 'COMMON.SETUP'
530 double precision emin
532 integer diff,allflag(maxdim),allnss,
533 & allihpb(maxdim),alljhpb(maxdim),
534 & newnss,newihpb(maxdim),newjhpb(maxdim)
536 integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
537 integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
542 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
551 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
555 if (allflag(i).eq.0 .and.
556 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
557 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
561 if (emin.lt.1.0d300) then
564 if (allflag(i).eq.0 .and.
565 & (allihpb(i).eq.allihpb(imin) .or.
566 & alljhpb(i).eq.allihpb(imin) .or.
567 & allihpb(i).eq.alljhpb(imin) .or.
568 & alljhpb(i).eq.alljhpb(imin))) then
575 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
579 if (allflag(i).eq.1) then
581 newihpb(newnss)=allihpb(i)
582 newjhpb(newnss)=alljhpb(i)
587 if (nfgtasks.gt.1)then
589 call MPI_Reduce(newnss,g_newnss,1,
590 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
591 call MPI_Gather(newnss,1,MPI_INTEGER,
592 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
595 displ(i)=i_newnss(i-1)+displ(i-1)
597 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
598 & g_newihpb,i_newnss,displ,MPI_INTEGER,
600 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
601 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
603 if(fg_rank.eq.0) then
604 c print *,'g_newnss',g_newnss
605 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
606 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
609 newihpb(i)=g_newihpb(i)
610 newjhpb(i)=g_newjhpb(i)
618 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
623 if (idssb(i).eq.newihpb(j) .and.
624 & jdssb(i).eq.newjhpb(j)) found=.true.
628 c if (.not.found.and.fg_rank.eq.0)
629 c & write(iout,'(a15,f12.2,f8.1,2i5)')
630 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
638 if (newihpb(i).eq.idssb(j) .and.
639 & newjhpb(i).eq.jdssb(j)) found=.true.
643 c if (.not.found.and.fg_rank.eq.0)
644 c & write(iout,'(a15,f12.2,f8.1,2i5)')
645 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
660 C-----------------------------------------------------------------------------
661 C-----------------------------------------------------------------------------
662 C-----------------------------------------------------------------------------
664 c$$$c-----------------------------------------------------------------------------
666 c$$$ subroutine ss_relax(i_in,j_in)
670 c$$$ include 'DIMENSIONS'
671 c$$$ include 'COMMON.VAR'
672 c$$$ include 'COMMON.CHAIN'
673 c$$$ include 'COMMON.IOUNITS'
674 c$$$ include 'COMMON.INTERACT'
676 c$$$c Input arguments
677 c$$$ integer i_in,j_in
679 c$$$c Local variables
680 c$$$ integer i,iretcode,nfun_sc
682 c$$$ double precision var(maxvar),e_sc,etot
689 c$$$ mask_side(i_in)=1
690 c$$$ mask_side(j_in)=1
692 c$$$c Minimize the two selected side-chains
693 c$$$ call overlap_sc(scfail) ! Better not fail!
694 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
701 c$$$c-------------------------------------------------------------
703 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
704 c$$$c Minimize side-chains only, starting from geom but without modifying
706 c$$$c If mask_r is already set, only the selected side-chains are minimized,
707 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
711 c$$$ include 'DIMENSIONS'
712 c$$$ include 'COMMON.IOUNITS'
713 c$$$ include 'COMMON.VAR'
714 c$$$ include 'COMMON.CHAIN'
715 c$$$ include 'COMMON.GEO'
716 c$$$ include 'COMMON.MINIM'
718 c$$$ common /srutu/ icall
720 c$$$c Output arguments
721 c$$$ double precision etot_sc
722 c$$$ integer iretcode,nfun
724 c$$$c External functions/subroutines
725 c$$$ external func_sc,grad_sc,fdum
727 c$$$c Local variables
729 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
731 c$$$ double precision rdum(1)
732 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
734 c$$$ integer i,nvar_restr
737 c$$$cmc start_minim=.true.
738 c$$$ call deflt(2,iv,liv,lv,v)
739 c$$$* 12 means fresh start, dont call deflt
741 c$$$* max num of fun calls
742 c$$$ if (maxfun.eq.0) maxfun=500
744 c$$$* max num of iterations
745 c$$$ if (maxmin.eq.0) maxmin=1000
747 c$$$* controls output
749 c$$$* selects output unit
751 c$$$c iv(21)=iout ! DEBUG
752 c$$$c iv(21)=8 ! DEBUG
753 c$$$* 1 means to print out result
755 c$$$c iv(22)=1 ! DEBUG
756 c$$$* 1 means to print out summary stats
758 c$$$c iv(23)=1 ! DEBUG
759 c$$$* 1 means to print initial x and d
761 c$$$c iv(24)=1 ! DEBUG
762 c$$$* min val for v(radfac) default is 0.1
764 c$$$* max val for v(radfac) default is 4.0
767 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
768 c$$$* the sumsl default is 0.1
770 c$$$* false conv if (act fnctn decrease) .lt. v(34)
771 c$$$* the sumsl default is 100*machep
772 c$$$ v(34)=v(34)/100.0D0
773 c$$$* absolute convergence
774 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
776 c$$$* relative convergence
777 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
779 c$$$* controls initial step size
781 c$$$* large vals of d correspond to small components of step
785 c$$$ do i=nphi+1,nvar
789 c$$$ call geom_to_var(nvar,x)
790 c$$$ IF (mask_r) THEN
791 c$$$ do i=1,nres ! Just in case...
795 c$$$ call x2xx(x,xx,nvar_restr)
796 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
797 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
800 c$$$c When minimizing ALL side-chains, etotal_sc is a little
801 c$$$c faster if we don't set mask_r
807 c$$$ call x2xx(x,xx,nvar_restr)
808 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
809 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
812 c$$$ call var_to_geom(nvar,x)
813 c$$$ call chainbuild_sc
820 c$$$C--------------------------------------------------------------------------
822 c$$$ subroutine chainbuild_sc
824 c$$$ include 'DIMENSIONS'
825 c$$$ include 'COMMON.VAR'
826 c$$$ include 'COMMON.INTERACT'
828 c$$$c Local variables
833 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
834 c$$$ call locate_side_chain(i)
841 c$$$C--------------------------------------------------------------------------
843 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
847 c$$$ include 'DIMENSIONS'
848 c$$$ include 'COMMON.DERIV'
849 c$$$ include 'COMMON.VAR'
850 c$$$ include 'COMMON.MINIM'
851 c$$$ include 'COMMON.IOUNITS'
853 c$$$c Input arguments
855 c$$$ double precision x(maxvar)
856 c$$$ double precision ufparm
859 c$$$c Input/Output arguments
861 c$$$ integer uiparm(1)
862 c$$$ double precision urparm(1)
864 c$$$c Output arguments
865 c$$$ double precision f
867 c$$$c Local variables
868 c$$$ double precision energia(0:n_ene)
870 c$$$c Variables used to intercept NaNs
871 c$$$ double precision x_sum
880 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
883 c$$$ x_sum=x_sum+x(i_NAN)
885 c$$$c Calculate the energy only if the coordinates are ok
886 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
887 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
893 c$$$ call var_to_geom_restr(n,x)
895 c$$$ call chainbuild_sc
896 c$$$ call etotal_sc(energia(0))
898 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
907 c$$$c-------------------------------------------------------
909 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
913 c$$$ include 'DIMENSIONS'
914 c$$$ include 'COMMON.CHAIN'
915 c$$$ include 'COMMON.DERIV'
916 c$$$ include 'COMMON.VAR'
917 c$$$ include 'COMMON.INTERACT'
918 c$$$ include 'COMMON.MINIM'
920 c$$$c Input arguments
922 c$$$ double precision x(maxvar)
923 c$$$ double precision ufparm
926 c$$$c Input/Output arguments
928 c$$$ integer uiparm(1)
929 c$$$ double precision urparm(1)
931 c$$$c Output arguments
932 c$$$ double precision g(maxvar)
934 c$$$c Local variables
935 c$$$ double precision f,gphii,gthetai,galphai,gomegai
936 c$$$ integer ig,ind,i,j,k,igall,ij
940 c$$$ if (nf-nfl+1) 20,30,40
941 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
942 c$$$c write (iout,*) 'grad 20'
943 c$$$ if (nf.eq.0) return
945 c$$$ 30 call var_to_geom_restr(n,x)
946 c$$$ call chainbuild_sc
948 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
952 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
958 c$$$ IF (mask_phi(i+2).eq.1) THEN
963 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
964 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
970 c$$$ ind=ind+nres-1-i
977 c$$$ IF (mask_theta(i+2).eq.1) THEN
983 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
984 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
989 c$$$ ind=ind+nres-1-i
994 c$$$ if (itype(i).ne.10) then
995 c$$$ IF (mask_side(i).eq.1) THEN
999 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1008 c$$$ if (itype(i).ne.10) then
1009 c$$$ IF (mask_side(i).eq.1) THEN
1013 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1021 c$$$C Add the components corresponding to local energy terms.
1028 c$$$ if (mask_phi(i).eq.1) then
1030 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1036 c$$$ if (mask_theta(i).eq.1) then
1038 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1044 c$$$ if (itype(i).ne.10) then
1046 c$$$ if (mask_side(i).eq.1) then
1048 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1055 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1061 c$$$C-----------------------------------------------------------------------------
1063 c$$$ subroutine etotal_sc(energy_sc)
1067 c$$$ include 'DIMENSIONS'
1068 c$$$ include 'COMMON.VAR'
1069 c$$$ include 'COMMON.INTERACT'
1070 c$$$ include 'COMMON.DERIV'
1071 c$$$ include 'COMMON.FFIELD'
1073 c$$$c Output arguments
1074 c$$$ double precision energy_sc(0:n_ene)
1076 c$$$c Local variables
1077 c$$$ double precision evdw,escloc
1082 c$$$ energy_sc(i)=0.0D0
1085 c$$$ if (mask_r) then
1086 c$$$ call egb_sc(evdw)
1087 c$$$ call esc_sc(escloc)
1090 c$$$ call esc(escloc)
1093 c$$$ if (evdw.eq.1.0D20) then
1094 c$$$ energy_sc(0)=evdw
1096 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1098 c$$$ energy_sc(1)=evdw
1099 c$$$ energy_sc(12)=escloc
1102 c$$$C Sum up the components of the Cartesian gradient.
1106 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1113 c$$$C-----------------------------------------------------------------------------
1115 c$$$ subroutine egb_sc(evdw)
1117 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1118 c$$$C assuming the Gay-Berne potential of interaction.
1120 c$$$ implicit real*8 (a-h,o-z)
1121 c$$$ include 'DIMENSIONS'
1122 c$$$ include 'COMMON.GEO'
1123 c$$$ include 'COMMON.VAR'
1124 c$$$ include 'COMMON.LOCAL'
1125 c$$$ include 'COMMON.CHAIN'
1126 c$$$ include 'COMMON.DERIV'
1127 c$$$ include 'COMMON.NAMES'
1128 c$$$ include 'COMMON.INTERACT'
1129 c$$$ include 'COMMON.IOUNITS'
1130 c$$$ include 'COMMON.CALC'
1131 c$$$ include 'COMMON.CONTROL'
1134 c$$$ energy_dec=.false.
1135 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1138 c$$$c if (icall.eq.0) lprn=.false.
1140 c$$$ do i=iatsc_s,iatsc_e
1142 c$$$ itypi1=itype(i+1)
1146 c$$$ dxi=dc_norm(1,nres+i)
1147 c$$$ dyi=dc_norm(2,nres+i)
1148 c$$$ dzi=dc_norm(3,nres+i)
1149 c$$$c dsci_inv=dsc_inv(itypi)
1150 c$$$ dsci_inv=vbld_inv(i+nres)
1151 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1152 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1154 c$$$C Calculate SC interaction energy.
1156 c$$$ do iint=1,nint_gr(i)
1157 c$$$ do j=istart(i,iint),iend(i,iint)
1158 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1161 c$$$c dscj_inv=dsc_inv(itypj)
1162 c$$$ dscj_inv=vbld_inv(j+nres)
1163 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1164 c$$$c & 1.0d0/vbld(j+nres)
1165 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1166 c$$$ sig0ij=sigma(itypi,itypj)
1167 c$$$ chi1=chi(itypi,itypj)
1168 c$$$ chi2=chi(itypj,itypi)
1169 c$$$ chi12=chi1*chi2
1170 c$$$ chip1=chip(itypi)
1171 c$$$ chip2=chip(itypj)
1172 c$$$ chip12=chip1*chip2
1173 c$$$ alf1=alp(itypi)
1174 c$$$ alf2=alp(itypj)
1175 c$$$ alf12=0.5D0*(alf1+alf2)
1176 c$$$C For diagnostics only!!!
1186 c$$$ xj=c(1,nres+j)-xi
1187 c$$$ yj=c(2,nres+j)-yi
1188 c$$$ zj=c(3,nres+j)-zi
1189 c$$$ dxj=dc_norm(1,nres+j)
1190 c$$$ dyj=dc_norm(2,nres+j)
1191 c$$$ dzj=dc_norm(3,nres+j)
1192 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1193 c$$$c write (iout,*) "j",j," dc_norm",
1194 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1195 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1196 c$$$ rij=dsqrt(rrij)
1197 c$$$C Calculate angle-dependent terms of energy and contributions to their
1199 c$$$ call sc_angular
1200 c$$$ sigsq=1.0D0/sigsq
1201 c$$$ sig=sig0ij*dsqrt(sigsq)
1202 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1203 c$$$c for diagnostics; uncomment
1204 c$$$c rij_shift=1.2*sig0ij
1205 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1206 c$$$ if (rij_shift.le.0.0D0) then
1208 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1209 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1210 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1213 c$$$ sigder=-sig*sigsq
1214 c$$$c---------------------------------------------------------------
1215 c$$$ rij_shift=1.0D0/rij_shift
1216 c$$$ fac=rij_shift**expon
1217 c$$$ e1=fac*fac*aa(itypi,itypj)
1218 c$$$ e2=fac*bb(itypi,itypj)
1219 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1220 c$$$ eps2der=evdwij*eps3rt
1221 c$$$ eps3der=evdwij*eps2rt
1222 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1223 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1224 c$$$ evdwij=evdwij*eps2rt*eps3rt
1225 c$$$ evdw=evdw+evdwij
1227 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1228 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1229 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1230 c$$$ & restyp(itypi),i,restyp(itypj),j,
1231 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1232 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1233 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1237 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1238 c$$$ & 'evdw',i,j,evdwij
1240 c$$$C Calculate gradient components.
1241 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1242 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1243 c$$$ sigder=fac*sigder
1246 c$$$C Calculate the radial part of the gradient
1250 c$$$C Calculate angular part of the gradient.
1256 c$$$ energy_dec=.false.
1260 c$$$c-----------------------------------------------------------------------------
1262 c$$$ subroutine esc_sc(escloc)
1263 c$$$C Calculate the local energy of a side chain and its derivatives in the
1264 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1265 c$$$C ALPHA and OMEGA.
1266 c$$$ implicit real*8 (a-h,o-z)
1267 c$$$ include 'DIMENSIONS'
1268 c$$$ include 'COMMON.GEO'
1269 c$$$ include 'COMMON.LOCAL'
1270 c$$$ include 'COMMON.VAR'
1271 c$$$ include 'COMMON.INTERACT'
1272 c$$$ include 'COMMON.DERIV'
1273 c$$$ include 'COMMON.CHAIN'
1274 c$$$ include 'COMMON.IOUNITS'
1275 c$$$ include 'COMMON.NAMES'
1276 c$$$ include 'COMMON.FFIELD'
1277 c$$$ include 'COMMON.CONTROL'
1278 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1279 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1280 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1281 c$$$ delta=0.02d0*pi
1283 c$$$c write (iout,'(a)') 'ESC'
1284 c$$$ do i=loc_start,loc_end
1285 c$$$ IF (mask_side(i).eq.1) THEN
1287 c$$$ if (it.eq.10) goto 1
1288 c$$$ nlobit=nlob(it)
1289 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1290 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1291 c$$$ theti=theta(i+1)-pipol
1292 c$$$ x(1)=dtan(theti)
1296 c$$$ if (x(2).gt.pi-delta) then
1298 c$$$ xtemp(2)=pi-delta
1300 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1302 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1303 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1304 c$$$ & escloci,dersc(2))
1305 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1306 c$$$ & ddersc0(1),dersc(1))
1307 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1308 c$$$ & ddersc0(3),dersc(3))
1309 c$$$ xtemp(2)=pi-delta
1310 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1312 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1313 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1314 c$$$ & dersc0(2),esclocbi,dersc02)
1315 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1316 c$$$ & dersc12,dersc01)
1317 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1318 c$$$ dersc0(1)=dersc01
1319 c$$$ dersc0(2)=dersc02
1320 c$$$ dersc0(3)=0.0d0
1322 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1324 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1325 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1326 c$$$c & esclocbi,ss,ssd
1327 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1328 c$$$c escloci=esclocbi
1329 c$$$c write (iout,*) escloci
1330 c$$$ else if (x(2).lt.delta) then
1334 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1336 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1337 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1338 c$$$ & escloci,dersc(2))
1339 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1340 c$$$ & ddersc0(1),dersc(1))
1341 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1342 c$$$ & ddersc0(3),dersc(3))
1344 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1346 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1347 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1348 c$$$ & dersc0(2),esclocbi,dersc02)
1349 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1350 c$$$ & dersc12,dersc01)
1351 c$$$ dersc0(1)=dersc01
1352 c$$$ dersc0(2)=dersc02
1353 c$$$ dersc0(3)=0.0d0
1354 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1356 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1358 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1359 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1360 c$$$c & esclocbi,ss,ssd
1361 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1362 c$$$c write (iout,*) escloci
1364 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1367 c$$$ escloc=escloc+escloci
1368 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1369 c$$$ & 'escloc',i,escloci
1370 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1372 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1373 c$$$ & wscloc*dersc(1)
1374 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1375 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1382 c$$$C-----------------------------------------------------------------------------
1384 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1386 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1387 c$$$C assuming the Gay-Berne potential of interaction.
1389 c$$$ implicit real*8 (a-h,o-z)
1390 c$$$ include 'DIMENSIONS'
1391 c$$$ include 'COMMON.GEO'
1392 c$$$ include 'COMMON.VAR'
1393 c$$$ include 'COMMON.LOCAL'
1394 c$$$ include 'COMMON.CHAIN'
1395 c$$$ include 'COMMON.DERIV'
1396 c$$$ include 'COMMON.NAMES'
1397 c$$$ include 'COMMON.INTERACT'
1398 c$$$ include 'COMMON.IOUNITS'
1399 c$$$ include 'COMMON.CALC'
1400 c$$$ include 'COMMON.CONTROL'
1403 c$$$ energy_dec=.false.
1404 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1408 c$$$c$$$ do i=iatsc_s,iatsc_e
1411 c$$$ itypi1=itype(i+1)
1415 c$$$ dxi=dc_norm(1,nres+i)
1416 c$$$ dyi=dc_norm(2,nres+i)
1417 c$$$ dzi=dc_norm(3,nres+i)
1418 c$$$c dsci_inv=dsc_inv(itypi)
1419 c$$$ dsci_inv=vbld_inv(i+nres)
1420 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1421 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1423 c$$$C Calculate SC interaction energy.
1425 c$$$c$$$ do iint=1,nint_gr(i)
1426 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1430 c$$$c dscj_inv=dsc_inv(itypj)
1431 c$$$ dscj_inv=vbld_inv(j+nres)
1432 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1433 c$$$c & 1.0d0/vbld(j+nres)
1434 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1435 c$$$ sig0ij=sigma(itypi,itypj)
1436 c$$$ chi1=chi(itypi,itypj)
1437 c$$$ chi2=chi(itypj,itypi)
1438 c$$$ chi12=chi1*chi2
1439 c$$$ chip1=chip(itypi)
1440 c$$$ chip2=chip(itypj)
1441 c$$$ chip12=chip1*chip2
1442 c$$$ alf1=alp(itypi)
1443 c$$$ alf2=alp(itypj)
1444 c$$$ alf12=0.5D0*(alf1+alf2)
1445 c$$$C For diagnostics only!!!
1455 c$$$ xj=c(1,nres+j)-xi
1456 c$$$ yj=c(2,nres+j)-yi
1457 c$$$ zj=c(3,nres+j)-zi
1458 c$$$ dxj=dc_norm(1,nres+j)
1459 c$$$ dyj=dc_norm(2,nres+j)
1460 c$$$ dzj=dc_norm(3,nres+j)
1461 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1462 c$$$c write (iout,*) "j",j," dc_norm",
1463 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1464 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1465 c$$$ rij=dsqrt(rrij)
1466 c$$$C Calculate angle-dependent terms of energy and contributions to their
1468 c$$$ call sc_angular
1469 c$$$ sigsq=1.0D0/sigsq
1470 c$$$ sig=sig0ij*dsqrt(sigsq)
1471 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1472 c$$$c for diagnostics; uncomment
1473 c$$$c rij_shift=1.2*sig0ij
1474 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1475 c$$$ if (rij_shift.le.0.0D0) then
1477 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1478 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1479 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1482 c$$$ sigder=-sig*sigsq
1483 c$$$c---------------------------------------------------------------
1484 c$$$ rij_shift=1.0D0/rij_shift
1485 c$$$ fac=rij_shift**expon
1486 c$$$ e1=fac*fac*aa(itypi,itypj)
1487 c$$$ e2=fac*bb(itypi,itypj)
1488 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1489 c$$$ eps2der=evdwij*eps3rt
1490 c$$$ eps3der=evdwij*eps2rt
1491 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1492 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1493 c$$$ evdwij=evdwij*eps2rt*eps3rt
1494 c$$$ evdw=evdw+evdwij
1496 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1497 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1498 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1499 c$$$ & restyp(itypi),i,restyp(itypj),j,
1500 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1501 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1502 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1506 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1507 c$$$ & 'evdw',i,j,evdwij
1509 c$$$C Calculate gradient components.
1510 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1511 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1512 c$$$ sigder=fac*sigder
1515 c$$$C Calculate the radial part of the gradient
1519 c$$$C Calculate angular part of the gradient.
1522 c$$$c$$$ enddo ! iint
1524 c$$$ energy_dec=.false.
1528 c$$$C-----------------------------------------------------------------------------
1530 c$$$ subroutine perturb_side_chain(i,angle)
1534 c$$$ include 'DIMENSIONS'
1535 c$$$ include 'COMMON.CHAIN'
1536 c$$$ include 'COMMON.GEO'
1537 c$$$ include 'COMMON.VAR'
1538 c$$$ include 'COMMON.LOCAL'
1539 c$$$ include 'COMMON.IOUNITS'
1541 c$$$c External functions
1542 c$$$ external ran_number
1543 c$$$ double precision ran_number
1545 c$$$c Input arguments
1547 c$$$ double precision angle ! In degrees
1549 c$$$c Local variables
1551 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1555 c$$$ rad_ang=angle*deg2rad
1558 c$$$ do while (length.lt.0.01)
1559 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1560 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1561 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1562 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1563 c$$$ + rand_v(3)*rand_v(3)
1564 c$$$ length=sqrt(length)
1565 c$$$ rand_v(1)=rand_v(1)/length
1566 c$$$ rand_v(2)=rand_v(2)/length
1567 c$$$ rand_v(3)=rand_v(3)/length
1568 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1569 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1570 c$$$ length=1.0D0-cost*cost
1571 c$$$ if (length.lt.0.0D0) length=0.0D0
1572 c$$$ length=sqrt(length)
1573 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1574 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1575 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1577 c$$$ rand_v(1)=rand_v(1)/length
1578 c$$$ rand_v(2)=rand_v(2)/length
1579 c$$$ rand_v(3)=rand_v(3)/length
1581 c$$$ cost=dcos(rad_ang)
1582 c$$$ sint=dsin(rad_ang)
1583 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1584 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1585 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1586 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1587 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1588 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1589 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1590 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1591 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1593 c$$$ call chainbuild_cart
1598 c$$$c----------------------------------------------------------------------------
1600 c$$$ subroutine ss_relax3(i_in,j_in)
1604 c$$$ include 'DIMENSIONS'
1605 c$$$ include 'COMMON.VAR'
1606 c$$$ include 'COMMON.CHAIN'
1607 c$$$ include 'COMMON.IOUNITS'
1608 c$$$ include 'COMMON.INTERACT'
1610 c$$$c External functions
1611 c$$$ external ran_number
1612 c$$$ double precision ran_number
1614 c$$$c Input arguments
1615 c$$$ integer i_in,j_in
1617 c$$$c Local variables
1618 c$$$ double precision energy_sc(0:n_ene),etot
1619 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1620 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1621 c$$$ integer n,i_pert,i
1622 c$$$ logical notdone
1631 c$$$ mask_side(i_in)=1
1632 c$$$ mask_side(j_in)=1
1634 c$$$ call etotal_sc(energy_sc)
1635 c$$$ etot=energy_sc(0)
1636 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1637 c$$$c + energy_sc(1),energy_sc(12)
1641 c$$$ do while (notdone)
1642 c$$$ if (mod(n,2).eq.0) then
1650 c$$$ org_dc(i)=dc(i,i_pert+nres)
1651 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1652 c$$$ org_c(i)=c(i,i_pert+nres)
1654 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1655 c$$$ call perturb_side_chain(i_pert,ang_pert)
1656 c$$$ call etotal_sc(energy_sc)
1657 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1658 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1659 c$$$ if (rand_fact.lt.exp_fact) then
1660 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1661 c$$$c + energy_sc(1),energy_sc(12)
1662 c$$$ etot=energy_sc(0)
1664 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1665 c$$$c + energy_sc(1),energy_sc(12)
1667 c$$$ dc(i,i_pert+nres)=org_dc(i)
1668 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1669 c$$$ c(i,i_pert+nres)=org_c(i)
1673 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1681 c$$$c----------------------------------------------------------------------------
1683 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1685 c$$$ include 'DIMENSIONS'
1687 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1688 c$$$*********************************************************************
1689 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1690 c$$$* the calling subprogram. *
1691 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1692 c$$$* calculated in the usual pythagorean way. *
1693 c$$$* absolute convergence occurs when the function is within v(31) of *
1694 c$$$* zero. unless you know the minimum value in advance, abs convg *
1695 c$$$* is probably not useful. *
1696 c$$$* relative convergence is when the model predicts that the function *
1697 c$$$* will decrease by less than v(32)*abs(fun). *
1698 c$$$*********************************************************************
1699 c$$$ include 'COMMON.IOUNITS'
1700 c$$$ include 'COMMON.VAR'
1701 c$$$ include 'COMMON.GEO'
1702 c$$$ include 'COMMON.MINIM'
1703 c$$$ include 'COMMON.CHAIN'
1705 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1706 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1707 c$$$ + orig_ss_dist(maxres2,maxres2)
1709 c$$$ double precision etot
1710 c$$$ integer iretcode,nfun,i_in,j_in
1713 c$$$ double precision dist
1714 c$$$ external ss_func,fdum
1715 c$$$ double precision ss_func,fdum
1717 c$$$ integer iv(liv),uiparm(2)
1718 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1722 c$$$ call deflt(2,iv,liv,lv,v)
1723 c$$$* 12 means fresh start, dont call deflt
1725 c$$$* max num of fun calls
1726 c$$$ if (maxfun.eq.0) maxfun=500
1728 c$$$* max num of iterations
1729 c$$$ if (maxmin.eq.0) maxmin=1000
1731 c$$$* controls output
1733 c$$$* selects output unit
1736 c$$$* 1 means to print out result
1738 c$$$* 1 means to print out summary stats
1740 c$$$* 1 means to print initial x and d
1742 c$$$* min val for v(radfac) default is 0.1
1744 c$$$* max val for v(radfac) default is 4.0
1747 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1748 c$$$* the sumsl default is 0.1
1750 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1751 c$$$* the sumsl default is 100*machep
1752 c$$$ v(34)=v(34)/100.0D0
1753 c$$$* absolute convergence
1754 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1757 c$$$* relative convergence
1758 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1761 c$$$* controls initial step size
1763 c$$$* large vals of d correspond to small components of step
1770 c$$$ orig_ss_dc(j,i)=dc(j,i)
1773 c$$$ call geom_to_var(nvar,orig_ss_var)
1777 c$$$ orig_ss_dist(j,i)=dist(j,i)
1778 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1779 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1780 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1792 c$$$ if (ialph(i,1).gt.0) then
1795 c$$$ x(k)=dc(j,i+nres)
1802 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1805 c$$$ nfun=iv(6)+iv(30)
1815 c$$$ if (ialph(i,1).gt.0) then
1818 c$$$ dc(j,i+nres)=x(k)
1822 c$$$ call chainbuild_cart
1827 c$$$C-----------------------------------------------------------------------------
1829 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1831 c$$$ include 'DIMENSIONS'
1832 c$$$ include 'COMMON.DERIV'
1833 c$$$ include 'COMMON.IOUNITS'
1834 c$$$ include 'COMMON.VAR'
1835 c$$$ include 'COMMON.CHAIN'
1836 c$$$ include 'COMMON.INTERACT'
1837 c$$$ include 'COMMON.SBRIDGE'
1839 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1840 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1841 c$$$ + orig_ss_dist(maxres2,maxres2)
1844 c$$$ double precision x(maxres6)
1846 c$$$ double precision f
1847 c$$$ integer uiparm(2)
1848 c$$$ real*8 urparm(1)
1849 c$$$ external ufparm
1850 c$$$ double precision ufparm
1853 c$$$ double precision dist
1855 c$$$ integer i,j,k,ss_i,ss_j
1856 c$$$ double precision tempf,var(maxvar)
1871 c$$$ if (ialph(i,1).gt.0) then
1874 c$$$ dc(j,i+nres)=x(k)
1878 c$$$ call chainbuild_cart
1880 c$$$ call geom_to_var(nvar,var)
1882 c$$$c Constraints on all angles
1884 c$$$ tempf=var(i)-orig_ss_var(i)
1885 c$$$ f=f+tempf*tempf
1888 c$$$c Constraints on all distances
1890 c$$$ if (i.gt.1) then
1891 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1892 c$$$ f=f+tempf*tempf
1895 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1896 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1897 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1898 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1899 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1900 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1901 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1902 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1906 c$$$c Constraints for the relevant CYS-CYS
1907 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1908 c$$$ f=f+tempf*tempf
1909 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
1911 c$$$c$$$ if (nf.ne.nfl) then
1912 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1913 c$$$c$$$ + f,dist(5+nres,14+nres)
1921 c$$$C-----------------------------------------------------------------------------
1922 c$$$C-----------------------------------------------------------------------------
1923 subroutine triple_ssbond_ene(resi,resj,resk,eij)
1924 include 'DIMENSIONS'
1925 include 'COMMON.SBRIDGE'
1926 include 'COMMON.CHAIN'
1927 include 'COMMON.DERIV'
1928 include 'COMMON.LOCAL'
1929 include 'COMMON.INTERACT'
1930 include 'COMMON.VAR'
1931 include 'COMMON.IOUNITS'
1932 include 'COMMON.CALC'
1939 c External functions
1940 double precision h_base
1944 integer resi,resj,resk
1947 double precision eij,eij1,eij2,eij3
1951 c integer itypi,itypj,k,l
1952 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
1953 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
1954 double precision xik,yik,zik,xjk,yjk,zjk
1955 double precision sig0ij,ljd,sig,fac,e1,e2
1956 double precision dcosom1(3),dcosom2(3),ed
1957 double precision pom1,pom2
1958 double precision ljA,ljB,ljXs
1959 double precision d_ljB(1:3)
1960 double precision ssA,ssB,ssC,ssXs
1961 double precision ssxm,ljxm,ssm,ljm
1962 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
1964 if (dtriss.eq.0) return
1968 C write(iout,*) resi,resj,resk
1970 dxi=dc_norm(1,nres+i)
1971 dyi=dc_norm(2,nres+i)
1972 dzi=dc_norm(3,nres+i)
1973 dsci_inv=vbld_inv(i+nres)
1983 dxj=dc_norm(1,nres+j)
1984 dyj=dc_norm(2,nres+j)
1985 dzj=dc_norm(3,nres+j)
1986 dscj_inv=vbld_inv(j+nres)
1992 dxk=dc_norm(1,nres+k)
1993 dyk=dc_norm(2,nres+k)
1994 dzk=dc_norm(3,nres+k)
1995 dscj_inv=vbld_inv(k+nres)
2005 rrij=(xij*xij+yij*yij+zij*zij)
2006 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2007 rrik=(xik*xik+yik*yik+zik*zik)
2009 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2011 C there are three combination of distances for each trisulfide bonds
2012 C The first case the ith atom is the center
2013 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2014 C distance y is second distance the a,b,c,d are parameters derived for
2015 C this problem d parameter was set as a penalty currenlty set to 1.
2016 if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then
2019 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss)
2021 C second case jth atom is center
2022 if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then
2025 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss)
2027 C the third case kth atom is the center
2028 if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then
2031 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss)
2037 C write(iout,*)i,j,k,eij
2038 C The energy penalty calculated now time for the gradient part
2039 C derivative over rij
2040 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5)
2041 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)
2046 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2047 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2050 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2051 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2053 C now derivative over rik
2054 fac=-eij1**2/dtriss*
2055 &(-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5)
2056 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
2061 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2062 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2065 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2066 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2068 C now derivative over rjk
2069 fac=-eij2**2/dtriss*
2070 &(-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)-
2071 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
2076 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2077 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2080 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2081 gvdwc(l,k)=gvdwc(l,k)+gg(l)