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(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)
256 eij=eij*eps2rt*eps3rt
259 e1=e1*eps1*eps2rt**2*eps3rt**2
260 ed=-expon*(e1+eij)/ljd
262 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
263 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
264 eom12=eij*eps1_om12+eps2der*eps2rt_om12
265 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
266 else if (rij.lt.ssxm) then
269 eij=ssA*ssd*ssd+ssB*ssd+ssC
271 ed=2*akcm*ssd+akct*deltat12
273 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
274 eom1=-2*akth*deltat1-pom1-om2*pom2
275 eom2= 2*akth*deltat2+pom1-om1*pom2
278 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
280 d_ssxm(1)=0.5D0*akct/ssA
284 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
285 d_ljxm(2)=d_ljxm(1)*sigsq_om2
286 d_ljxm(3)=d_ljxm(1)*sigsq_om12
287 d_ljxm(1)=d_ljxm(1)*sigsq_om1
289 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
292 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
296 ssm=ssC-0.25D0*ssB*ssB/ssA
297 d_ssm(1)=0.5D0*akct*ssB/ssA
298 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
299 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
301 f1=(rij-xm)/(ssxm-xm)
302 f2=(rij-ssxm)/(xm-ssxm)
306 delta_inv=1.0d0/(xm-ssxm)
307 deltasq_inv=delta_inv*delta_inv
309 fac1=deltasq_inv*fac*(xm-rij)
310 fac2=deltasq_inv*fac*(rij-ssxm)
311 ed=delta_inv*(Ht*hd2-ssm*hd1)
312 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
313 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
314 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
317 ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
318 d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
319 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
320 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
322 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
323 f1=(rij-ljxm)/(xm-ljxm)
324 f2=(rij-xm)/(ljxm-xm)
328 delta_inv=1.0d0/(ljxm-xm)
329 deltasq_inv=delta_inv*delta_inv
331 fac1=deltasq_inv*fac*(ljxm-rij)
332 fac2=deltasq_inv*fac*(rij-xm)
333 ed=delta_inv*(ljm*hd2-Ht*hd1)
334 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
335 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
336 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
338 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
340 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
346 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
347 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
348 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
350 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
351 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
352 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
353 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
356 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
358 c$$$ d_ljm(k)=ljm*d_ljB(k)
362 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
363 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
364 c$$$ d_ss(2)=akct*ssd
365 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
366 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
369 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
370 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
371 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
373 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
374 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
376 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
378 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
379 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
380 c$$$ h1=h_base(f1,hd1)
381 c$$$ h2=h_base(f2,hd2)
382 c$$$ eij=ss*h1+ljf*h2
383 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
384 c$$$ deltasq_inv=delta_inv*delta_inv
385 c$$$ fac=ljf*hd2-ss*hd1
386 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
387 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
388 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
389 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
390 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
391 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
392 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
394 c$$$ havebond=.false.
395 c$$$ if (ed.gt.0.0d0) havebond=.true.
396 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
403 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
404 c write(iout,'(a15,f12.2,f8.1,2i5)')
405 c & "SSBOND_E_FORM",totT,t_bath,i,j
409 dyn_ssbond_ij(i,j)=eij
410 else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
411 dyn_ssbond_ij(i,j)=1.0d300
414 c write(iout,'(a15,f12.2,f8.1,2i5)')
415 c & "SSBOND_E_BREAK",totT,t_bath,i,j
422 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
423 & "CHECKSTOP",rij,eij,ed
428 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
435 c-------END TESTING CODE
438 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
439 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
442 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
445 gvdwx(k,i)=gvdwx(k,i)-gg(k)
446 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
447 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
448 gvdwx(k,j)=gvdwx(k,j)+gg(k)
449 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
450 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
454 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
459 gvdwc(l,i)=gvdwc(l,i)-gg(l)
460 gvdwc(l,j)=gvdwc(l,j)+gg(l)
466 C-----------------------------------------------------------------------------
468 double precision function h_base(x,deriv)
469 c A smooth function going 0->1 in range [0,1]
470 c It should NOT be called outside range [0,1], it will not work there.
477 double precision deriv
483 c Two parabolas put together. First derivative zero at extrema
484 c$$$ if (x.lt.0.5D0) then
485 c$$$ h_base=2.0D0*x*x
489 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
490 c$$$ deriv=4.0D0*deriv
493 c Third degree polynomial. First derivative zero at extrema
494 h_base=x*x*(3.0d0-2.0d0*x)
495 deriv=6.0d0*x*(1.0d0-x)
497 c Fifth degree polynomial. First and second derivatives zero at extrema
499 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
501 c$$$ deriv=deriv*deriv
502 c$$$ deriv=30.0d0*xsq*deriv
507 c----------------------------------------------------------------------------
509 subroutine dyn_set_nss
510 c Adjust nss and other relevant variables based on dyn_ssbond_ij
518 include 'COMMON.SBRIDGE'
519 include 'COMMON.CHAIN'
520 include 'COMMON.IOUNITS'
521 include 'COMMON.SETUP'
529 double precision emin
531 integer diff,allflag(maxdim),allnss,
532 & allihpb(maxdim),alljhpb(maxdim),
533 & newnss,newihpb(maxdim),newjhpb(maxdim)
535 integer i_newnss(max_fg_procs),displ(max_fg_procs)
536 integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
541 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
550 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
554 if (allflag(i).eq.0 .and.
555 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
556 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
560 if (emin.lt.1.0d300) then
563 if (allflag(i).eq.0 .and.
564 & (allihpb(i).eq.allihpb(imin) .or.
565 & alljhpb(i).eq.allihpb(imin) .or.
566 & allihpb(i).eq.alljhpb(imin) .or.
567 & alljhpb(i).eq.alljhpb(imin))) then
574 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
578 if (allflag(i).eq.1) then
580 newihpb(newnss)=allihpb(i)
581 newjhpb(newnss)=alljhpb(i)
586 if (nfgtasks.gt.1)then
588 call MPI_Reduce(newnss,g_newnss,1,
589 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
590 call MPI_Gather(newnss,1,MPI_INTEGER,
591 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
594 displ(i)=i_newnss(i-1)+displ(i-1)
596 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
597 & g_newihpb,i_newnss,displ,MPI_INTEGER,
599 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
600 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
602 if(fg_rank.eq.0) then
603 c print *,'g_newnss',g_newnss
604 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
605 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
608 newihpb(i)=g_newihpb(i)
609 newjhpb(i)=g_newjhpb(i)
617 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
622 if (idssb(i).eq.newihpb(j) .and.
623 & jdssb(i).eq.newjhpb(j)) found=.true.
627 if (.not.found.and.fg_rank.eq.0)
628 & write(iout,'(a15,f12.2,f8.1,2i5)')
629 & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
637 if (newihpb(i).eq.idssb(j) .and.
638 & newjhpb(i).eq.jdssb(j)) found=.true.
642 if (.not.found.and.fg_rank.eq.0)
643 & write(iout,'(a15,f12.2,f8.1,2i5)')
644 & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
658 c----------------------------------------------------------------------------
661 subroutine read_ssHist
666 include "DIMENSIONS.FREE"
667 include 'COMMON.FREE'
671 character*80 controlcard
674 call card_concat(controlcard,.true.)
676 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
683 c----------------------------------------------------------------------------
686 C-----------------------------------------------------------------------------
687 C-----------------------------------------------------------------------------
688 C-----------------------------------------------------------------------------
689 C-----------------------------------------------------------------------------
690 C-----------------------------------------------------------------------------
691 C-----------------------------------------------------------------------------
692 C-----------------------------------------------------------------------------
694 c$$$c-----------------------------------------------------------------------------
696 c$$$ subroutine ss_relax(i_in,j_in)
700 c$$$ include 'DIMENSIONS'
701 c$$$ include 'COMMON.VAR'
702 c$$$ include 'COMMON.CHAIN'
703 c$$$ include 'COMMON.IOUNITS'
704 c$$$ include 'COMMON.INTERACT'
706 c$$$c Input arguments
707 c$$$ integer i_in,j_in
709 c$$$c Local variables
710 c$$$ integer i,iretcode,nfun_sc
712 c$$$ double precision var(maxvar),e_sc,etot
719 c$$$ mask_side(i_in)=1
720 c$$$ mask_side(j_in)=1
722 c$$$c Minimize the two selected side-chains
723 c$$$ call overlap_sc(scfail) ! Better not fail!
724 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
731 c$$$c-------------------------------------------------------------
733 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
734 c$$$c Minimize side-chains only, starting from geom but without modifying
736 c$$$c If mask_r is already set, only the selected side-chains are minimized,
737 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
741 c$$$ include 'DIMENSIONS'
742 c$$$ include 'COMMON.IOUNITS'
743 c$$$ include 'COMMON.VAR'
744 c$$$ include 'COMMON.CHAIN'
745 c$$$ include 'COMMON.GEO'
746 c$$$ include 'COMMON.MINIM'
748 c$$$ common /srutu/ icall
750 c$$$c Output arguments
751 c$$$ double precision etot_sc
752 c$$$ integer iretcode,nfun
754 c$$$c External functions/subroutines
755 c$$$ external func_sc,grad_sc,fdum
757 c$$$c Local variables
759 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
761 c$$$ double precision rdum(1)
762 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
764 c$$$ integer i,nvar_restr
767 c$$$cmc start_minim=.true.
768 c$$$ call deflt(2,iv,liv,lv,v)
769 c$$$* 12 means fresh start, dont call deflt
771 c$$$* max num of fun calls
772 c$$$ if (maxfun.eq.0) maxfun=500
774 c$$$* max num of iterations
775 c$$$ if (maxmin.eq.0) maxmin=1000
777 c$$$* controls output
779 c$$$* selects output unit
781 c$$$c iv(21)=iout ! DEBUG
782 c$$$c iv(21)=8 ! DEBUG
783 c$$$* 1 means to print out result
785 c$$$c iv(22)=1 ! DEBUG
786 c$$$* 1 means to print out summary stats
788 c$$$c iv(23)=1 ! DEBUG
789 c$$$* 1 means to print initial x and d
791 c$$$c iv(24)=1 ! DEBUG
792 c$$$* min val for v(radfac) default is 0.1
794 c$$$* max val for v(radfac) default is 4.0
797 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
798 c$$$* the sumsl default is 0.1
800 c$$$* false conv if (act fnctn decrease) .lt. v(34)
801 c$$$* the sumsl default is 100*machep
802 c$$$ v(34)=v(34)/100.0D0
803 c$$$* absolute convergence
804 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
806 c$$$* relative convergence
807 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
809 c$$$* controls initial step size
811 c$$$* large vals of d correspond to small components of step
815 c$$$ do i=nphi+1,nvar
819 c$$$ call geom_to_var(nvar,x)
820 c$$$ IF (mask_r) THEN
821 c$$$ do i=1,nres ! Just in case...
825 c$$$ call x2xx(x,xx,nvar_restr)
826 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
827 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
830 c$$$c When minimizing ALL side-chains, etotal_sc is a little
831 c$$$c faster if we don't set mask_r
837 c$$$ call x2xx(x,xx,nvar_restr)
838 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
839 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
842 c$$$ call var_to_geom(nvar,x)
843 c$$$ call chainbuild_sc
850 c$$$C--------------------------------------------------------------------------
852 c$$$ subroutine chainbuild_sc
854 c$$$ include 'DIMENSIONS'
855 c$$$ include 'COMMON.VAR'
856 c$$$ include 'COMMON.INTERACT'
858 c$$$c Local variables
863 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
864 c$$$ call locate_side_chain(i)
871 c$$$C--------------------------------------------------------------------------
873 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
877 c$$$ include 'DIMENSIONS'
878 c$$$ include 'COMMON.DERIV'
879 c$$$ include 'COMMON.VAR'
880 c$$$ include 'COMMON.MINIM'
881 c$$$ include 'COMMON.IOUNITS'
883 c$$$c Input arguments
885 c$$$ double precision x(maxvar)
886 c$$$ double precision ufparm
889 c$$$c Input/Output arguments
891 c$$$ integer uiparm(1)
892 c$$$ double precision urparm(1)
894 c$$$c Output arguments
895 c$$$ double precision f
897 c$$$c Local variables
898 c$$$ double precision energia(0:n_ene)
900 c$$$c Variables used to intercept NaNs
901 c$$$ double precision x_sum
910 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
913 c$$$ x_sum=x_sum+x(i_NAN)
915 c$$$c Calculate the energy only if the coordinates are ok
916 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
917 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
923 c$$$ call var_to_geom_restr(n,x)
925 c$$$ call chainbuild_sc
926 c$$$ call etotal_sc(energia(0))
928 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
937 c$$$c-------------------------------------------------------
939 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
943 c$$$ include 'DIMENSIONS'
944 c$$$ include 'COMMON.CHAIN'
945 c$$$ include 'COMMON.DERIV'
946 c$$$ include 'COMMON.VAR'
947 c$$$ include 'COMMON.INTERACT'
948 c$$$ include 'COMMON.MINIM'
950 c$$$c Input arguments
952 c$$$ double precision x(maxvar)
953 c$$$ double precision ufparm
956 c$$$c Input/Output arguments
958 c$$$ integer uiparm(1)
959 c$$$ double precision urparm(1)
961 c$$$c Output arguments
962 c$$$ double precision g(maxvar)
964 c$$$c Local variables
965 c$$$ double precision f,gphii,gthetai,galphai,gomegai
966 c$$$ integer ig,ind,i,j,k,igall,ij
970 c$$$ if (nf-nfl+1) 20,30,40
971 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
972 c$$$c write (iout,*) 'grad 20'
973 c$$$ if (nf.eq.0) return
975 c$$$ 30 call var_to_geom_restr(n,x)
976 c$$$ call chainbuild_sc
978 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
982 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
988 c$$$ IF (mask_phi(i+2).eq.1) THEN
993 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
994 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1000 c$$$ ind=ind+nres-1-i
1007 c$$$ IF (mask_theta(i+2).eq.1) THEN
1010 c$$$ do j=i+1,nres-1
1013 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1014 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1019 c$$$ ind=ind+nres-1-i
1024 c$$$ if (itype(i).ne.10) then
1025 c$$$ IF (mask_side(i).eq.1) THEN
1029 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1038 c$$$ if (itype(i).ne.10) then
1039 c$$$ IF (mask_side(i).eq.1) THEN
1043 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1051 c$$$C Add the components corresponding to local energy terms.
1058 c$$$ if (mask_phi(i).eq.1) then
1060 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1066 c$$$ if (mask_theta(i).eq.1) then
1068 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1074 c$$$ if (itype(i).ne.10) then
1076 c$$$ if (mask_side(i).eq.1) then
1078 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1085 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1091 c$$$C-----------------------------------------------------------------------------
1093 c$$$ subroutine etotal_sc(energy_sc)
1097 c$$$ include 'DIMENSIONS'
1098 c$$$ include 'COMMON.VAR'
1099 c$$$ include 'COMMON.INTERACT'
1100 c$$$ include 'COMMON.DERIV'
1101 c$$$ include 'COMMON.FFIELD'
1103 c$$$c Output arguments
1104 c$$$ double precision energy_sc(0:n_ene)
1106 c$$$c Local variables
1107 c$$$ double precision evdw,escloc
1112 c$$$ energy_sc(i)=0.0D0
1115 c$$$ if (mask_r) then
1116 c$$$ call egb_sc(evdw)
1117 c$$$ call esc_sc(escloc)
1120 c$$$ call esc(escloc)
1123 c$$$ if (evdw.eq.1.0D20) then
1124 c$$$ energy_sc(0)=evdw
1126 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1128 c$$$ energy_sc(1)=evdw
1129 c$$$ energy_sc(12)=escloc
1132 c$$$C Sum up the components of the Cartesian gradient.
1136 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1143 c$$$C-----------------------------------------------------------------------------
1145 c$$$ subroutine egb_sc(evdw)
1147 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1148 c$$$C assuming the Gay-Berne potential of interaction.
1150 c$$$ implicit real*8 (a-h,o-z)
1151 c$$$ include 'DIMENSIONS'
1152 c$$$ include 'COMMON.GEO'
1153 c$$$ include 'COMMON.VAR'
1154 c$$$ include 'COMMON.LOCAL'
1155 c$$$ include 'COMMON.CHAIN'
1156 c$$$ include 'COMMON.DERIV'
1157 c$$$ include 'COMMON.NAMES'
1158 c$$$ include 'COMMON.INTERACT'
1159 c$$$ include 'COMMON.IOUNITS'
1160 c$$$ include 'COMMON.CALC'
1161 c$$$ include 'COMMON.CONTROL'
1164 c$$$ energy_dec=.false.
1165 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1168 c$$$c if (icall.eq.0) lprn=.false.
1170 c$$$ do i=iatsc_s,iatsc_e
1172 c$$$ itypi1=itype(i+1)
1176 c$$$ dxi=dc_norm(1,nres+i)
1177 c$$$ dyi=dc_norm(2,nres+i)
1178 c$$$ dzi=dc_norm(3,nres+i)
1179 c$$$c dsci_inv=dsc_inv(itypi)
1180 c$$$ dsci_inv=vbld_inv(i+nres)
1181 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1182 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1184 c$$$C Calculate SC interaction energy.
1186 c$$$ do iint=1,nint_gr(i)
1187 c$$$ do j=istart(i,iint),iend(i,iint)
1188 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1191 c$$$c dscj_inv=dsc_inv(itypj)
1192 c$$$ dscj_inv=vbld_inv(j+nres)
1193 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1194 c$$$c & 1.0d0/vbld(j+nres)
1195 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1196 c$$$ sig0ij=sigma(itypi,itypj)
1197 c$$$ chi1=chi(itypi,itypj)
1198 c$$$ chi2=chi(itypj,itypi)
1199 c$$$ chi12=chi1*chi2
1200 c$$$ chip1=chip(itypi)
1201 c$$$ chip2=chip(itypj)
1202 c$$$ chip12=chip1*chip2
1203 c$$$ alf1=alp(itypi)
1204 c$$$ alf2=alp(itypj)
1205 c$$$ alf12=0.5D0*(alf1+alf2)
1206 c$$$C For diagnostics only!!!
1216 c$$$ xj=c(1,nres+j)-xi
1217 c$$$ yj=c(2,nres+j)-yi
1218 c$$$ zj=c(3,nres+j)-zi
1219 c$$$ dxj=dc_norm(1,nres+j)
1220 c$$$ dyj=dc_norm(2,nres+j)
1221 c$$$ dzj=dc_norm(3,nres+j)
1222 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1223 c$$$c write (iout,*) "j",j," dc_norm",
1224 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1225 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1226 c$$$ rij=dsqrt(rrij)
1227 c$$$C Calculate angle-dependent terms of energy and contributions to their
1229 c$$$ call sc_angular
1230 c$$$ sigsq=1.0D0/sigsq
1231 c$$$ sig=sig0ij*dsqrt(sigsq)
1232 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1233 c$$$c for diagnostics; uncomment
1234 c$$$c rij_shift=1.2*sig0ij
1235 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1236 c$$$ if (rij_shift.le.0.0D0) then
1238 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1239 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1240 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1243 c$$$ sigder=-sig*sigsq
1244 c$$$c---------------------------------------------------------------
1245 c$$$ rij_shift=1.0D0/rij_shift
1246 c$$$ fac=rij_shift**expon
1247 c$$$ e1=fac*fac*aa(itypi,itypj)
1248 c$$$ e2=fac*bb(itypi,itypj)
1249 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1250 c$$$ eps2der=evdwij*eps3rt
1251 c$$$ eps3der=evdwij*eps2rt
1252 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1253 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1254 c$$$ evdwij=evdwij*eps2rt*eps3rt
1255 c$$$ evdw=evdw+evdwij
1257 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1258 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1259 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1260 c$$$ & restyp(itypi),i,restyp(itypj),j,
1261 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1262 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1263 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1267 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1268 c$$$ & 'evdw',i,j,evdwij
1270 c$$$C Calculate gradient components.
1271 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1272 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1273 c$$$ sigder=fac*sigder
1276 c$$$C Calculate the radial part of the gradient
1280 c$$$C Calculate angular part of the gradient.
1286 c$$$ energy_dec=.false.
1290 c$$$c-----------------------------------------------------------------------------
1292 c$$$ subroutine esc_sc(escloc)
1293 c$$$C Calculate the local energy of a side chain and its derivatives in the
1294 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1295 c$$$C ALPHA and OMEGA.
1296 c$$$ implicit real*8 (a-h,o-z)
1297 c$$$ include 'DIMENSIONS'
1298 c$$$ include 'COMMON.GEO'
1299 c$$$ include 'COMMON.LOCAL'
1300 c$$$ include 'COMMON.VAR'
1301 c$$$ include 'COMMON.INTERACT'
1302 c$$$ include 'COMMON.DERIV'
1303 c$$$ include 'COMMON.CHAIN'
1304 c$$$ include 'COMMON.IOUNITS'
1305 c$$$ include 'COMMON.NAMES'
1306 c$$$ include 'COMMON.FFIELD'
1307 c$$$ include 'COMMON.CONTROL'
1308 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1309 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1310 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1311 c$$$ delta=0.02d0*pi
1313 c$$$c write (iout,'(a)') 'ESC'
1314 c$$$ do i=loc_start,loc_end
1315 c$$$ IF (mask_side(i).eq.1) THEN
1317 c$$$ if (it.eq.10) goto 1
1318 c$$$ nlobit=nlob(it)
1319 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1320 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1321 c$$$ theti=theta(i+1)-pipol
1322 c$$$ x(1)=dtan(theti)
1326 c$$$ if (x(2).gt.pi-delta) then
1328 c$$$ xtemp(2)=pi-delta
1330 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1332 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1333 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1334 c$$$ & escloci,dersc(2))
1335 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1336 c$$$ & ddersc0(1),dersc(1))
1337 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1338 c$$$ & ddersc0(3),dersc(3))
1339 c$$$ xtemp(2)=pi-delta
1340 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1342 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1343 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1344 c$$$ & dersc0(2),esclocbi,dersc02)
1345 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1346 c$$$ & dersc12,dersc01)
1347 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1348 c$$$ dersc0(1)=dersc01
1349 c$$$ dersc0(2)=dersc02
1350 c$$$ dersc0(3)=0.0d0
1352 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1354 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1355 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1356 c$$$c & esclocbi,ss,ssd
1357 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1358 c$$$c escloci=esclocbi
1359 c$$$c write (iout,*) escloci
1360 c$$$ else if (x(2).lt.delta) then
1364 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1366 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1367 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1368 c$$$ & escloci,dersc(2))
1369 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1370 c$$$ & ddersc0(1),dersc(1))
1371 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1372 c$$$ & ddersc0(3),dersc(3))
1374 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1376 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1377 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1378 c$$$ & dersc0(2),esclocbi,dersc02)
1379 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1380 c$$$ & dersc12,dersc01)
1381 c$$$ dersc0(1)=dersc01
1382 c$$$ dersc0(2)=dersc02
1383 c$$$ dersc0(3)=0.0d0
1384 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1386 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1388 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1389 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1390 c$$$c & esclocbi,ss,ssd
1391 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1392 c$$$c write (iout,*) escloci
1394 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1397 c$$$ escloc=escloc+escloci
1398 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1399 c$$$ & 'escloc',i,escloci
1400 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1402 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1403 c$$$ & wscloc*dersc(1)
1404 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1405 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1412 c$$$C-----------------------------------------------------------------------------
1414 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1416 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1417 c$$$C assuming the Gay-Berne potential of interaction.
1419 c$$$ implicit real*8 (a-h,o-z)
1420 c$$$ include 'DIMENSIONS'
1421 c$$$ include 'COMMON.GEO'
1422 c$$$ include 'COMMON.VAR'
1423 c$$$ include 'COMMON.LOCAL'
1424 c$$$ include 'COMMON.CHAIN'
1425 c$$$ include 'COMMON.DERIV'
1426 c$$$ include 'COMMON.NAMES'
1427 c$$$ include 'COMMON.INTERACT'
1428 c$$$ include 'COMMON.IOUNITS'
1429 c$$$ include 'COMMON.CALC'
1430 c$$$ include 'COMMON.CONTROL'
1433 c$$$ energy_dec=.false.
1434 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1438 c$$$c$$$ do i=iatsc_s,iatsc_e
1441 c$$$ itypi1=itype(i+1)
1445 c$$$ dxi=dc_norm(1,nres+i)
1446 c$$$ dyi=dc_norm(2,nres+i)
1447 c$$$ dzi=dc_norm(3,nres+i)
1448 c$$$c dsci_inv=dsc_inv(itypi)
1449 c$$$ dsci_inv=vbld_inv(i+nres)
1450 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1451 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1453 c$$$C Calculate SC interaction energy.
1455 c$$$c$$$ do iint=1,nint_gr(i)
1456 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1460 c$$$c dscj_inv=dsc_inv(itypj)
1461 c$$$ dscj_inv=vbld_inv(j+nres)
1462 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1463 c$$$c & 1.0d0/vbld(j+nres)
1464 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1465 c$$$ sig0ij=sigma(itypi,itypj)
1466 c$$$ chi1=chi(itypi,itypj)
1467 c$$$ chi2=chi(itypj,itypi)
1468 c$$$ chi12=chi1*chi2
1469 c$$$ chip1=chip(itypi)
1470 c$$$ chip2=chip(itypj)
1471 c$$$ chip12=chip1*chip2
1472 c$$$ alf1=alp(itypi)
1473 c$$$ alf2=alp(itypj)
1474 c$$$ alf12=0.5D0*(alf1+alf2)
1475 c$$$C For diagnostics only!!!
1485 c$$$ xj=c(1,nres+j)-xi
1486 c$$$ yj=c(2,nres+j)-yi
1487 c$$$ zj=c(3,nres+j)-zi
1488 c$$$ dxj=dc_norm(1,nres+j)
1489 c$$$ dyj=dc_norm(2,nres+j)
1490 c$$$ dzj=dc_norm(3,nres+j)
1491 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1492 c$$$c write (iout,*) "j",j," dc_norm",
1493 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1494 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1495 c$$$ rij=dsqrt(rrij)
1496 c$$$C Calculate angle-dependent terms of energy and contributions to their
1498 c$$$ call sc_angular
1499 c$$$ sigsq=1.0D0/sigsq
1500 c$$$ sig=sig0ij*dsqrt(sigsq)
1501 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1502 c$$$c for diagnostics; uncomment
1503 c$$$c rij_shift=1.2*sig0ij
1504 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1505 c$$$ if (rij_shift.le.0.0D0) then
1507 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1508 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1509 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1512 c$$$ sigder=-sig*sigsq
1513 c$$$c---------------------------------------------------------------
1514 c$$$ rij_shift=1.0D0/rij_shift
1515 c$$$ fac=rij_shift**expon
1516 c$$$ e1=fac*fac*aa(itypi,itypj)
1517 c$$$ e2=fac*bb(itypi,itypj)
1518 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1519 c$$$ eps2der=evdwij*eps3rt
1520 c$$$ eps3der=evdwij*eps2rt
1521 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1522 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1523 c$$$ evdwij=evdwij*eps2rt*eps3rt
1524 c$$$ evdw=evdw+evdwij
1526 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1527 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1528 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1529 c$$$ & restyp(itypi),i,restyp(itypj),j,
1530 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1531 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1532 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1536 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1537 c$$$ & 'evdw',i,j,evdwij
1539 c$$$C Calculate gradient components.
1540 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1541 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1542 c$$$ sigder=fac*sigder
1545 c$$$C Calculate the radial part of the gradient
1549 c$$$C Calculate angular part of the gradient.
1552 c$$$c$$$ enddo ! iint
1554 c$$$ energy_dec=.false.
1558 c$$$C-----------------------------------------------------------------------------
1560 c$$$ subroutine perturb_side_chain(i,angle)
1564 c$$$ include 'DIMENSIONS'
1565 c$$$ include 'COMMON.CHAIN'
1566 c$$$ include 'COMMON.GEO'
1567 c$$$ include 'COMMON.VAR'
1568 c$$$ include 'COMMON.LOCAL'
1569 c$$$ include 'COMMON.IOUNITS'
1571 c$$$c External functions
1572 c$$$ external ran_number
1573 c$$$ double precision ran_number
1575 c$$$c Input arguments
1577 c$$$ double precision angle ! In degrees
1579 c$$$c Local variables
1581 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1585 c$$$ rad_ang=angle*deg2rad
1588 c$$$ do while (length.lt.0.01)
1589 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1590 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1591 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1592 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1593 c$$$ + rand_v(3)*rand_v(3)
1594 c$$$ length=sqrt(length)
1595 c$$$ rand_v(1)=rand_v(1)/length
1596 c$$$ rand_v(2)=rand_v(2)/length
1597 c$$$ rand_v(3)=rand_v(3)/length
1598 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1599 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1600 c$$$ length=1.0D0-cost*cost
1601 c$$$ if (length.lt.0.0D0) length=0.0D0
1602 c$$$ length=sqrt(length)
1603 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1604 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1605 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1607 c$$$ rand_v(1)=rand_v(1)/length
1608 c$$$ rand_v(2)=rand_v(2)/length
1609 c$$$ rand_v(3)=rand_v(3)/length
1611 c$$$ cost=dcos(rad_ang)
1612 c$$$ sint=dsin(rad_ang)
1613 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1614 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1615 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1616 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1617 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1618 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1619 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1620 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1621 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1623 c$$$ call chainbuild_cart
1628 c$$$c----------------------------------------------------------------------------
1630 c$$$ subroutine ss_relax3(i_in,j_in)
1634 c$$$ include 'DIMENSIONS'
1635 c$$$ include 'COMMON.VAR'
1636 c$$$ include 'COMMON.CHAIN'
1637 c$$$ include 'COMMON.IOUNITS'
1638 c$$$ include 'COMMON.INTERACT'
1640 c$$$c External functions
1641 c$$$ external ran_number
1642 c$$$ double precision ran_number
1644 c$$$c Input arguments
1645 c$$$ integer i_in,j_in
1647 c$$$c Local variables
1648 c$$$ double precision energy_sc(0:n_ene),etot
1649 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1650 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1651 c$$$ integer n,i_pert,i
1652 c$$$ logical notdone
1661 c$$$ mask_side(i_in)=1
1662 c$$$ mask_side(j_in)=1
1664 c$$$ call etotal_sc(energy_sc)
1665 c$$$ etot=energy_sc(0)
1666 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1667 c$$$c + energy_sc(1),energy_sc(12)
1671 c$$$ do while (notdone)
1672 c$$$ if (mod(n,2).eq.0) then
1680 c$$$ org_dc(i)=dc(i,i_pert+nres)
1681 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1682 c$$$ org_c(i)=c(i,i_pert+nres)
1684 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1685 c$$$ call perturb_side_chain(i_pert,ang_pert)
1686 c$$$ call etotal_sc(energy_sc)
1687 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1688 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1689 c$$$ if (rand_fact.lt.exp_fact) then
1690 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1691 c$$$c + energy_sc(1),energy_sc(12)
1692 c$$$ etot=energy_sc(0)
1694 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1695 c$$$c + energy_sc(1),energy_sc(12)
1697 c$$$ dc(i,i_pert+nres)=org_dc(i)
1698 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1699 c$$$ c(i,i_pert+nres)=org_c(i)
1703 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1711 c$$$c----------------------------------------------------------------------------
1713 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1715 c$$$ include 'DIMENSIONS'
1717 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1718 c$$$*********************************************************************
1719 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1720 c$$$* the calling subprogram. *
1721 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1722 c$$$* calculated in the usual pythagorean way. *
1723 c$$$* absolute convergence occurs when the function is within v(31) of *
1724 c$$$* zero. unless you know the minimum value in advance, abs convg *
1725 c$$$* is probably not useful. *
1726 c$$$* relative convergence is when the model predicts that the function *
1727 c$$$* will decrease by less than v(32)*abs(fun). *
1728 c$$$*********************************************************************
1729 c$$$ include 'COMMON.IOUNITS'
1730 c$$$ include 'COMMON.VAR'
1731 c$$$ include 'COMMON.GEO'
1732 c$$$ include 'COMMON.MINIM'
1733 c$$$ include 'COMMON.CHAIN'
1735 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1736 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1737 c$$$ + orig_ss_dist(maxres2,maxres2)
1739 c$$$ double precision etot
1740 c$$$ integer iretcode,nfun,i_in,j_in
1743 c$$$ double precision dist
1744 c$$$ external ss_func,fdum
1745 c$$$ double precision ss_func,fdum
1747 c$$$ integer iv(liv),uiparm(2)
1748 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1752 c$$$ call deflt(2,iv,liv,lv,v)
1753 c$$$* 12 means fresh start, dont call deflt
1755 c$$$* max num of fun calls
1756 c$$$ if (maxfun.eq.0) maxfun=500
1758 c$$$* max num of iterations
1759 c$$$ if (maxmin.eq.0) maxmin=1000
1761 c$$$* controls output
1763 c$$$* selects output unit
1766 c$$$* 1 means to print out result
1768 c$$$* 1 means to print out summary stats
1770 c$$$* 1 means to print initial x and d
1772 c$$$* min val for v(radfac) default is 0.1
1774 c$$$* max val for v(radfac) default is 4.0
1777 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1778 c$$$* the sumsl default is 0.1
1780 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1781 c$$$* the sumsl default is 100*machep
1782 c$$$ v(34)=v(34)/100.0D0
1783 c$$$* absolute convergence
1784 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1787 c$$$* relative convergence
1788 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1791 c$$$* controls initial step size
1793 c$$$* large vals of d correspond to small components of step
1800 c$$$ orig_ss_dc(j,i)=dc(j,i)
1803 c$$$ call geom_to_var(nvar,orig_ss_var)
1807 c$$$ orig_ss_dist(j,i)=dist(j,i)
1808 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1809 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1810 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1822 c$$$ if (ialph(i,1).gt.0) then
1825 c$$$ x(k)=dc(j,i+nres)
1832 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1835 c$$$ nfun=iv(6)+iv(30)
1845 c$$$ if (ialph(i,1).gt.0) then
1848 c$$$ dc(j,i+nres)=x(k)
1852 c$$$ call chainbuild_cart
1857 c$$$C-----------------------------------------------------------------------------
1859 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1861 c$$$ include 'DIMENSIONS'
1862 c$$$ include 'COMMON.DERIV'
1863 c$$$ include 'COMMON.IOUNITS'
1864 c$$$ include 'COMMON.VAR'
1865 c$$$ include 'COMMON.CHAIN'
1866 c$$$ include 'COMMON.INTERACT'
1867 c$$$ include 'COMMON.SBRIDGE'
1869 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1870 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1871 c$$$ + orig_ss_dist(maxres2,maxres2)
1874 c$$$ double precision x(maxres6)
1876 c$$$ double precision f
1877 c$$$ integer uiparm(2)
1878 c$$$ real*8 urparm(1)
1879 c$$$ external ufparm
1880 c$$$ double precision ufparm
1883 c$$$ double precision dist
1885 c$$$ integer i,j,k,ss_i,ss_j
1886 c$$$ double precision tempf,var(maxvar)
1901 c$$$ if (ialph(i,1).gt.0) then
1904 c$$$ dc(j,i+nres)=x(k)
1908 c$$$ call chainbuild_cart
1910 c$$$ call geom_to_var(nvar,var)
1912 c$$$c Constraints on all angles
1914 c$$$ tempf=var(i)-orig_ss_var(i)
1915 c$$$ f=f+tempf*tempf
1918 c$$$c Constraints on all distances
1920 c$$$ if (i.gt.1) then
1921 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1922 c$$$ f=f+tempf*tempf
1925 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1926 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1927 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1928 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1929 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1930 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1931 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1932 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1936 c$$$c Constraints for the relevant CYS-CYS
1937 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1938 c$$$ f=f+tempf*tempf
1939 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
1941 c$$$c$$$ if (nf.ne.nfl) then
1942 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1943 c$$$c$$$ + f,dist(5+nres,14+nres)
1951 c$$$C-----------------------------------------------------------------------------
1952 subroutine triple_ssbond_ene(resi,resj,resk,eij)
1953 include 'DIMENSIONS'
1954 include 'COMMON.SBRIDGE'
1955 include 'COMMON.CHAIN'
1956 include 'COMMON.DERIV'
1957 include 'COMMON.LOCAL'
1958 include 'COMMON.INTERACT'
1959 include 'COMMON.VAR'
1960 include 'COMMON.IOUNITS'
1961 include 'COMMON.CALC'
1968 c External functions
1969 double precision h_base
1973 integer resi,resj,resk
1976 double precision eij,eij1,eij2,eij3
1980 c integer itypi,itypj,k,l
1981 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
1982 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
1983 double precision xik,yik,zik,xjk,yjk,zjk
1984 double precision sig0ij,ljd,sig,fac,e1,e2
1985 double precision dcosom1(3),dcosom2(3),ed
1986 double precision pom1,pom2
1987 double precision ljA,ljB,ljXs
1988 double precision d_ljB(1:3)
1989 double precision ssA,ssB,ssC,ssXs
1990 double precision ssxm,ljxm,ssm,ljm
1991 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
1996 C write(iout,*) resi,resj,resk
1998 dxi=dc_norm(1,nres+i)
1999 dyi=dc_norm(2,nres+i)
2000 dzi=dc_norm(3,nres+i)
2001 dsci_inv=vbld_inv(i+nres)
2011 dxj=dc_norm(1,nres+j)
2012 dyj=dc_norm(2,nres+j)
2013 dzj=dc_norm(3,nres+j)
2014 dscj_inv=vbld_inv(j+nres)
2020 dxk=dc_norm(1,nres+k)
2021 dyk=dc_norm(2,nres+k)
2022 dzk=dc_norm(3,nres+k)
2023 dscj_inv=vbld_inv(k+nres)
2033 rrij=(xij*xij+yij*yij+zij*zij)
2034 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2035 rrik=(xik*xik+yik*yik+zik*zik)
2037 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2039 C there are three combination of distances for each trisulfide bonds
2040 C The first case the ith atom is the center
2041 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2042 C distance y is second distance the a,b,c,d are parameters derived for
2043 C this problem d parameter was set as a penalty currenlty set to 1.
2044 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2045 C second case jth atom is center
2046 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2047 C the third case kth atom is the center
2048 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2053 C write(iout,*)i,j,k,eij
2054 C The energy penalty calculated now time for the gradient part
2055 C derivative over rij
2056 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2057 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2062 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2063 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2066 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2067 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2069 C now derivative over rik
2070 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2071 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2076 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2077 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2080 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2081 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2083 C now derivative over rjk
2084 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2085 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2090 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2091 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2094 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2095 gvdwc(l,k)=gvdwc(l,k)+gg(l)