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 C 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)
662 c$$$c-----------------------------------------------------------------------------
664 c$$$ subroutine ss_relax(i_in,j_in)
668 c$$$ include 'DIMENSIONS'
669 c$$$ include 'COMMON.VAR'
670 c$$$ include 'COMMON.CHAIN'
671 c$$$ include 'COMMON.IOUNITS'
672 c$$$ include 'COMMON.INTERACT'
674 c$$$c Input arguments
675 c$$$ integer i_in,j_in
677 c$$$c Local variables
678 c$$$ integer i,iretcode,nfun_sc
680 c$$$ double precision var(maxvar),e_sc,etot
687 c$$$ mask_side(i_in)=1
688 c$$$ mask_side(j_in)=1
690 c$$$c Minimize the two selected side-chains
691 c$$$ call overlap_sc(scfail) ! Better not fail!
692 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
699 c$$$c-------------------------------------------------------------
701 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
702 c$$$c Minimize side-chains only, starting from geom but without modifying
704 c$$$c If mask_r is already set, only the selected side-chains are minimized,
705 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
709 c$$$ include 'DIMENSIONS'
710 c$$$ include 'COMMON.IOUNITS'
711 c$$$ include 'COMMON.VAR'
712 c$$$ include 'COMMON.CHAIN'
713 c$$$ include 'COMMON.GEO'
714 c$$$ include 'COMMON.MINIM'
716 c$$$ common /srutu/ icall
718 c$$$c Output arguments
719 c$$$ double precision etot_sc
720 c$$$ integer iretcode,nfun
722 c$$$c External functions/subroutines
723 c$$$ external func_sc,grad_sc,fdum
725 c$$$c Local variables
727 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
729 c$$$ double precision rdum(1)
730 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
732 c$$$ integer i,nvar_restr
735 c$$$cmc start_minim=.true.
736 c$$$ call deflt(2,iv,liv,lv,v)
737 c$$$* 12 means fresh start, dont call deflt
739 c$$$* max num of fun calls
740 c$$$ if (maxfun.eq.0) maxfun=500
742 c$$$* max num of iterations
743 c$$$ if (maxmin.eq.0) maxmin=1000
745 c$$$* controls output
747 c$$$* selects output unit
749 c$$$c iv(21)=iout ! DEBUG
750 c$$$c iv(21)=8 ! DEBUG
751 c$$$* 1 means to print out result
753 c$$$c iv(22)=1 ! DEBUG
754 c$$$* 1 means to print out summary stats
756 c$$$c iv(23)=1 ! DEBUG
757 c$$$* 1 means to print initial x and d
759 c$$$c iv(24)=1 ! DEBUG
760 c$$$* min val for v(radfac) default is 0.1
762 c$$$* max val for v(radfac) default is 4.0
765 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
766 c$$$* the sumsl default is 0.1
768 c$$$* false conv if (act fnctn decrease) .lt. v(34)
769 c$$$* the sumsl default is 100*machep
770 c$$$ v(34)=v(34)/100.0D0
771 c$$$* absolute convergence
772 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
774 c$$$* relative convergence
775 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
777 c$$$* controls initial step size
779 c$$$* large vals of d correspond to small components of step
783 c$$$ do i=nphi+1,nvar
787 c$$$ call geom_to_var(nvar,x)
788 c$$$ IF (mask_r) THEN
789 c$$$ do i=1,nres ! Just in case...
793 c$$$ call x2xx(x,xx,nvar_restr)
794 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
795 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
798 c$$$c When minimizing ALL side-chains, etotal_sc is a little
799 c$$$c faster if we don't set mask_r
805 c$$$ call x2xx(x,xx,nvar_restr)
806 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
807 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
810 c$$$ call var_to_geom(nvar,x)
811 c$$$ call chainbuild_sc
818 c$$$C--------------------------------------------------------------------------
820 c$$$ subroutine chainbuild_sc
822 c$$$ include 'DIMENSIONS'
823 c$$$ include 'COMMON.VAR'
824 c$$$ include 'COMMON.INTERACT'
826 c$$$c Local variables
831 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
832 c$$$ call locate_side_chain(i)
839 c$$$C--------------------------------------------------------------------------
841 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
845 c$$$ include 'DIMENSIONS'
846 c$$$ include 'COMMON.DERIV'
847 c$$$ include 'COMMON.VAR'
848 c$$$ include 'COMMON.MINIM'
849 c$$$ include 'COMMON.IOUNITS'
851 c$$$c Input arguments
853 c$$$ double precision x(maxvar)
854 c$$$ double precision ufparm
857 c$$$c Input/Output arguments
859 c$$$ integer uiparm(1)
860 c$$$ double precision urparm(1)
862 c$$$c Output arguments
863 c$$$ double precision f
865 c$$$c Local variables
866 c$$$ double precision energia(0:n_ene)
868 c$$$c Variables used to intercept NaNs
869 c$$$ double precision x_sum
878 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
881 c$$$ x_sum=x_sum+x(i_NAN)
883 c$$$c Calculate the energy only if the coordinates are ok
884 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
885 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
891 c$$$ call var_to_geom_restr(n,x)
893 c$$$ call chainbuild_sc
894 c$$$ call etotal_sc(energia(0))
896 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
905 c$$$c-------------------------------------------------------
907 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
911 c$$$ include 'DIMENSIONS'
912 c$$$ include 'COMMON.CHAIN'
913 c$$$ include 'COMMON.DERIV'
914 c$$$ include 'COMMON.VAR'
915 c$$$ include 'COMMON.INTERACT'
916 c$$$ include 'COMMON.MINIM'
918 c$$$c Input arguments
920 c$$$ double precision x(maxvar)
921 c$$$ double precision ufparm
924 c$$$c Input/Output arguments
926 c$$$ integer uiparm(1)
927 c$$$ double precision urparm(1)
929 c$$$c Output arguments
930 c$$$ double precision g(maxvar)
932 c$$$c Local variables
933 c$$$ double precision f,gphii,gthetai,galphai,gomegai
934 c$$$ integer ig,ind,i,j,k,igall,ij
938 c$$$ if (nf-nfl+1) 20,30,40
939 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
940 c$$$c write (iout,*) 'grad 20'
941 c$$$ if (nf.eq.0) return
943 c$$$ 30 call var_to_geom_restr(n,x)
944 c$$$ call chainbuild_sc
946 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
950 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
956 c$$$ IF (mask_phi(i+2).eq.1) THEN
961 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
962 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
968 c$$$ ind=ind+nres-1-i
975 c$$$ IF (mask_theta(i+2).eq.1) THEN
981 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
982 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
987 c$$$ ind=ind+nres-1-i
992 c$$$ if (itype(i).ne.10) then
993 c$$$ IF (mask_side(i).eq.1) THEN
997 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1006 c$$$ if (itype(i).ne.10) then
1007 c$$$ IF (mask_side(i).eq.1) THEN
1011 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1019 c$$$C Add the components corresponding to local energy terms.
1026 c$$$ if (mask_phi(i).eq.1) then
1028 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1034 c$$$ if (mask_theta(i).eq.1) then
1036 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1042 c$$$ if (itype(i).ne.10) then
1044 c$$$ if (mask_side(i).eq.1) then
1046 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1053 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1059 c$$$C-----------------------------------------------------------------------------
1061 c$$$ subroutine etotal_sc(energy_sc)
1065 c$$$ include 'DIMENSIONS'
1066 c$$$ include 'COMMON.VAR'
1067 c$$$ include 'COMMON.INTERACT'
1068 c$$$ include 'COMMON.DERIV'
1069 c$$$ include 'COMMON.FFIELD'
1071 c$$$c Output arguments
1072 c$$$ double precision energy_sc(0:n_ene)
1074 c$$$c Local variables
1075 c$$$ double precision evdw,escloc
1080 c$$$ energy_sc(i)=0.0D0
1083 c$$$ if (mask_r) then
1084 c$$$ call egb_sc(evdw)
1085 c$$$ call esc_sc(escloc)
1088 c$$$ call esc(escloc)
1091 c$$$ if (evdw.eq.1.0D20) then
1092 c$$$ energy_sc(0)=evdw
1094 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1096 c$$$ energy_sc(1)=evdw
1097 c$$$ energy_sc(12)=escloc
1100 c$$$C Sum up the components of the Cartesian gradient.
1104 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1111 c$$$C-----------------------------------------------------------------------------
1113 c$$$ subroutine egb_sc(evdw)
1115 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1116 c$$$C assuming the Gay-Berne potential of interaction.
1118 c$$$ implicit real*8 (a-h,o-z)
1119 c$$$ include 'DIMENSIONS'
1120 c$$$ include 'COMMON.GEO'
1121 c$$$ include 'COMMON.VAR'
1122 c$$$ include 'COMMON.LOCAL'
1123 c$$$ include 'COMMON.CHAIN'
1124 c$$$ include 'COMMON.DERIV'
1125 c$$$ include 'COMMON.NAMES'
1126 c$$$ include 'COMMON.INTERACT'
1127 c$$$ include 'COMMON.IOUNITS'
1128 c$$$ include 'COMMON.CALC'
1129 c$$$ include 'COMMON.CONTROL'
1132 c$$$ energy_dec=.false.
1133 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1136 c$$$c if (icall.eq.0) lprn=.false.
1138 c$$$ do i=iatsc_s,iatsc_e
1140 c$$$ itypi1=itype(i+1)
1144 c$$$ dxi=dc_norm(1,nres+i)
1145 c$$$ dyi=dc_norm(2,nres+i)
1146 c$$$ dzi=dc_norm(3,nres+i)
1147 c$$$c dsci_inv=dsc_inv(itypi)
1148 c$$$ dsci_inv=vbld_inv(i+nres)
1149 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1150 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1152 c$$$C Calculate SC interaction energy.
1154 c$$$ do iint=1,nint_gr(i)
1155 c$$$ do j=istart(i,iint),iend(i,iint)
1156 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1159 c$$$c dscj_inv=dsc_inv(itypj)
1160 c$$$ dscj_inv=vbld_inv(j+nres)
1161 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1162 c$$$c & 1.0d0/vbld(j+nres)
1163 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1164 c$$$ sig0ij=sigma(itypi,itypj)
1165 c$$$ chi1=chi(itypi,itypj)
1166 c$$$ chi2=chi(itypj,itypi)
1167 c$$$ chi12=chi1*chi2
1168 c$$$ chip1=chip(itypi)
1169 c$$$ chip2=chip(itypj)
1170 c$$$ chip12=chip1*chip2
1171 c$$$ alf1=alp(itypi)
1172 c$$$ alf2=alp(itypj)
1173 c$$$ alf12=0.5D0*(alf1+alf2)
1174 c$$$C For diagnostics only!!!
1184 c$$$ xj=c(1,nres+j)-xi
1185 c$$$ yj=c(2,nres+j)-yi
1186 c$$$ zj=c(3,nres+j)-zi
1187 c$$$ dxj=dc_norm(1,nres+j)
1188 c$$$ dyj=dc_norm(2,nres+j)
1189 c$$$ dzj=dc_norm(3,nres+j)
1190 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1191 c$$$c write (iout,*) "j",j," dc_norm",
1192 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1193 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1194 c$$$ rij=dsqrt(rrij)
1195 c$$$C Calculate angle-dependent terms of energy and contributions to their
1197 c$$$ call sc_angular
1198 c$$$ sigsq=1.0D0/sigsq
1199 c$$$ sig=sig0ij*dsqrt(sigsq)
1200 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1201 c$$$c for diagnostics; uncomment
1202 c$$$c rij_shift=1.2*sig0ij
1203 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1204 c$$$ if (rij_shift.le.0.0D0) then
1206 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1207 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1208 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1211 c$$$ sigder=-sig*sigsq
1212 c$$$c---------------------------------------------------------------
1213 c$$$ rij_shift=1.0D0/rij_shift
1214 c$$$ fac=rij_shift**expon
1215 c$$$ e1=fac*fac*aa(itypi,itypj)
1216 c$$$ e2=fac*bb(itypi,itypj)
1217 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1218 c$$$ eps2der=evdwij*eps3rt
1219 c$$$ eps3der=evdwij*eps2rt
1220 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1221 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1222 c$$$ evdwij=evdwij*eps2rt*eps3rt
1223 c$$$ evdw=evdw+evdwij
1225 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1226 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1227 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1228 c$$$ & restyp(itypi),i,restyp(itypj),j,
1229 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1230 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1231 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1235 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1236 c$$$ & 'evdw',i,j,evdwij
1238 c$$$C Calculate gradient components.
1239 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1240 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1241 c$$$ sigder=fac*sigder
1244 c$$$C Calculate the radial part of the gradient
1248 c$$$C Calculate angular part of the gradient.
1254 c$$$ energy_dec=.false.
1258 c$$$c-----------------------------------------------------------------------------
1260 c$$$ subroutine esc_sc(escloc)
1261 c$$$C Calculate the local energy of a side chain and its derivatives in the
1262 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1263 c$$$C ALPHA and OMEGA.
1264 c$$$ implicit real*8 (a-h,o-z)
1265 c$$$ include 'DIMENSIONS'
1266 c$$$ include 'COMMON.GEO'
1267 c$$$ include 'COMMON.LOCAL'
1268 c$$$ include 'COMMON.VAR'
1269 c$$$ include 'COMMON.INTERACT'
1270 c$$$ include 'COMMON.DERIV'
1271 c$$$ include 'COMMON.CHAIN'
1272 c$$$ include 'COMMON.IOUNITS'
1273 c$$$ include 'COMMON.NAMES'
1274 c$$$ include 'COMMON.FFIELD'
1275 c$$$ include 'COMMON.CONTROL'
1276 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1277 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1278 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1279 c$$$ delta=0.02d0*pi
1281 c$$$c write (iout,'(a)') 'ESC'
1282 c$$$ do i=loc_start,loc_end
1283 c$$$ IF (mask_side(i).eq.1) THEN
1285 c$$$ if (it.eq.10) goto 1
1286 c$$$ nlobit=nlob(it)
1287 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1288 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1289 c$$$ theti=theta(i+1)-pipol
1290 c$$$ x(1)=dtan(theti)
1294 c$$$ if (x(2).gt.pi-delta) then
1296 c$$$ xtemp(2)=pi-delta
1298 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1300 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1301 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1302 c$$$ & escloci,dersc(2))
1303 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1304 c$$$ & ddersc0(1),dersc(1))
1305 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1306 c$$$ & ddersc0(3),dersc(3))
1307 c$$$ xtemp(2)=pi-delta
1308 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1310 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1311 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1312 c$$$ & dersc0(2),esclocbi,dersc02)
1313 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1314 c$$$ & dersc12,dersc01)
1315 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1316 c$$$ dersc0(1)=dersc01
1317 c$$$ dersc0(2)=dersc02
1318 c$$$ dersc0(3)=0.0d0
1320 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1322 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1323 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1324 c$$$c & esclocbi,ss,ssd
1325 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1326 c$$$c escloci=esclocbi
1327 c$$$c write (iout,*) escloci
1328 c$$$ else if (x(2).lt.delta) then
1332 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1334 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1335 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1336 c$$$ & escloci,dersc(2))
1337 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1338 c$$$ & ddersc0(1),dersc(1))
1339 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1340 c$$$ & ddersc0(3),dersc(3))
1342 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1344 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1345 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1346 c$$$ & dersc0(2),esclocbi,dersc02)
1347 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1348 c$$$ & dersc12,dersc01)
1349 c$$$ dersc0(1)=dersc01
1350 c$$$ dersc0(2)=dersc02
1351 c$$$ dersc0(3)=0.0d0
1352 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1354 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1356 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1357 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1358 c$$$c & esclocbi,ss,ssd
1359 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1360 c$$$c write (iout,*) escloci
1362 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1365 c$$$ escloc=escloc+escloci
1366 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1367 c$$$ & 'escloc',i,escloci
1368 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1370 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1371 c$$$ & wscloc*dersc(1)
1372 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1373 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1380 c$$$C-----------------------------------------------------------------------------
1382 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1384 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1385 c$$$C assuming the Gay-Berne potential of interaction.
1387 c$$$ implicit real*8 (a-h,o-z)
1388 c$$$ include 'DIMENSIONS'
1389 c$$$ include 'COMMON.GEO'
1390 c$$$ include 'COMMON.VAR'
1391 c$$$ include 'COMMON.LOCAL'
1392 c$$$ include 'COMMON.CHAIN'
1393 c$$$ include 'COMMON.DERIV'
1394 c$$$ include 'COMMON.NAMES'
1395 c$$$ include 'COMMON.INTERACT'
1396 c$$$ include 'COMMON.IOUNITS'
1397 c$$$ include 'COMMON.CALC'
1398 c$$$ include 'COMMON.CONTROL'
1401 c$$$ energy_dec=.false.
1402 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1406 c$$$c$$$ do i=iatsc_s,iatsc_e
1409 c$$$ itypi1=itype(i+1)
1413 c$$$ dxi=dc_norm(1,nres+i)
1414 c$$$ dyi=dc_norm(2,nres+i)
1415 c$$$ dzi=dc_norm(3,nres+i)
1416 c$$$c dsci_inv=dsc_inv(itypi)
1417 c$$$ dsci_inv=vbld_inv(i+nres)
1418 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1419 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1421 c$$$C Calculate SC interaction energy.
1423 c$$$c$$$ do iint=1,nint_gr(i)
1424 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1428 c$$$c dscj_inv=dsc_inv(itypj)
1429 c$$$ dscj_inv=vbld_inv(j+nres)
1430 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1431 c$$$c & 1.0d0/vbld(j+nres)
1432 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1433 c$$$ sig0ij=sigma(itypi,itypj)
1434 c$$$ chi1=chi(itypi,itypj)
1435 c$$$ chi2=chi(itypj,itypi)
1436 c$$$ chi12=chi1*chi2
1437 c$$$ chip1=chip(itypi)
1438 c$$$ chip2=chip(itypj)
1439 c$$$ chip12=chip1*chip2
1440 c$$$ alf1=alp(itypi)
1441 c$$$ alf2=alp(itypj)
1442 c$$$ alf12=0.5D0*(alf1+alf2)
1443 c$$$C For diagnostics only!!!
1453 c$$$ xj=c(1,nres+j)-xi
1454 c$$$ yj=c(2,nres+j)-yi
1455 c$$$ zj=c(3,nres+j)-zi
1456 c$$$ dxj=dc_norm(1,nres+j)
1457 c$$$ dyj=dc_norm(2,nres+j)
1458 c$$$ dzj=dc_norm(3,nres+j)
1459 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1460 c$$$c write (iout,*) "j",j," dc_norm",
1461 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1462 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1463 c$$$ rij=dsqrt(rrij)
1464 c$$$C Calculate angle-dependent terms of energy and contributions to their
1466 c$$$ call sc_angular
1467 c$$$ sigsq=1.0D0/sigsq
1468 c$$$ sig=sig0ij*dsqrt(sigsq)
1469 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1470 c$$$c for diagnostics; uncomment
1471 c$$$c rij_shift=1.2*sig0ij
1472 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1473 c$$$ if (rij_shift.le.0.0D0) then
1475 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1476 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1477 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1480 c$$$ sigder=-sig*sigsq
1481 c$$$c---------------------------------------------------------------
1482 c$$$ rij_shift=1.0D0/rij_shift
1483 c$$$ fac=rij_shift**expon
1484 c$$$ e1=fac*fac*aa(itypi,itypj)
1485 c$$$ e2=fac*bb(itypi,itypj)
1486 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1487 c$$$ eps2der=evdwij*eps3rt
1488 c$$$ eps3der=evdwij*eps2rt
1489 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1490 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1491 c$$$ evdwij=evdwij*eps2rt*eps3rt
1492 c$$$ evdw=evdw+evdwij
1494 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1495 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1496 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1497 c$$$ & restyp(itypi),i,restyp(itypj),j,
1498 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1499 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1500 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1504 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1505 c$$$ & 'evdw',i,j,evdwij
1507 c$$$C Calculate gradient components.
1508 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1509 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1510 c$$$ sigder=fac*sigder
1513 c$$$C Calculate the radial part of the gradient
1517 c$$$C Calculate angular part of the gradient.
1520 c$$$c$$$ enddo ! iint
1522 c$$$ energy_dec=.false.
1526 c$$$C-----------------------------------------------------------------------------
1528 c$$$ subroutine perturb_side_chain(i,angle)
1532 c$$$ include 'DIMENSIONS'
1533 c$$$ include 'COMMON.CHAIN'
1534 c$$$ include 'COMMON.GEO'
1535 c$$$ include 'COMMON.VAR'
1536 c$$$ include 'COMMON.LOCAL'
1537 c$$$ include 'COMMON.IOUNITS'
1539 c$$$c External functions
1540 c$$$ external ran_number
1541 c$$$ double precision ran_number
1543 c$$$c Input arguments
1545 c$$$ double precision angle ! In degrees
1547 c$$$c Local variables
1549 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1553 c$$$ rad_ang=angle*deg2rad
1556 c$$$ do while (length.lt.0.01)
1557 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1558 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1559 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1560 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1561 c$$$ + rand_v(3)*rand_v(3)
1562 c$$$ length=sqrt(length)
1563 c$$$ rand_v(1)=rand_v(1)/length
1564 c$$$ rand_v(2)=rand_v(2)/length
1565 c$$$ rand_v(3)=rand_v(3)/length
1566 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1567 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1568 c$$$ length=1.0D0-cost*cost
1569 c$$$ if (length.lt.0.0D0) length=0.0D0
1570 c$$$ length=sqrt(length)
1571 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1572 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1573 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1575 c$$$ rand_v(1)=rand_v(1)/length
1576 c$$$ rand_v(2)=rand_v(2)/length
1577 c$$$ rand_v(3)=rand_v(3)/length
1579 c$$$ cost=dcos(rad_ang)
1580 c$$$ sint=dsin(rad_ang)
1581 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1582 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1583 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1584 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1585 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1586 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1587 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1588 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1589 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1591 c$$$ call chainbuild_cart
1596 c$$$c----------------------------------------------------------------------------
1598 c$$$ subroutine ss_relax3(i_in,j_in)
1602 c$$$ include 'DIMENSIONS'
1603 c$$$ include 'COMMON.VAR'
1604 c$$$ include 'COMMON.CHAIN'
1605 c$$$ include 'COMMON.IOUNITS'
1606 c$$$ include 'COMMON.INTERACT'
1608 c$$$c External functions
1609 c$$$ external ran_number
1610 c$$$ double precision ran_number
1612 c$$$c Input arguments
1613 c$$$ integer i_in,j_in
1615 c$$$c Local variables
1616 c$$$ double precision energy_sc(0:n_ene),etot
1617 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1618 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1619 c$$$ integer n,i_pert,i
1620 c$$$ logical notdone
1629 c$$$ mask_side(i_in)=1
1630 c$$$ mask_side(j_in)=1
1632 c$$$ call etotal_sc(energy_sc)
1633 c$$$ etot=energy_sc(0)
1634 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1635 c$$$c + energy_sc(1),energy_sc(12)
1639 c$$$ do while (notdone)
1640 c$$$ if (mod(n,2).eq.0) then
1648 c$$$ org_dc(i)=dc(i,i_pert+nres)
1649 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1650 c$$$ org_c(i)=c(i,i_pert+nres)
1652 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1653 c$$$ call perturb_side_chain(i_pert,ang_pert)
1654 c$$$ call etotal_sc(energy_sc)
1655 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1656 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1657 c$$$ if (rand_fact.lt.exp_fact) then
1658 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1659 c$$$c + energy_sc(1),energy_sc(12)
1660 c$$$ etot=energy_sc(0)
1662 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1663 c$$$c + energy_sc(1),energy_sc(12)
1665 c$$$ dc(i,i_pert+nres)=org_dc(i)
1666 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1667 c$$$ c(i,i_pert+nres)=org_c(i)
1671 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1679 c$$$c----------------------------------------------------------------------------
1681 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1683 c$$$ include 'DIMENSIONS'
1685 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1686 c$$$*********************************************************************
1687 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1688 c$$$* the calling subprogram. *
1689 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1690 c$$$* calculated in the usual pythagorean way. *
1691 c$$$* absolute convergence occurs when the function is within v(31) of *
1692 c$$$* zero. unless you know the minimum value in advance, abs convg *
1693 c$$$* is probably not useful. *
1694 c$$$* relative convergence is when the model predicts that the function *
1695 c$$$* will decrease by less than v(32)*abs(fun). *
1696 c$$$*********************************************************************
1697 c$$$ include 'COMMON.IOUNITS'
1698 c$$$ include 'COMMON.VAR'
1699 c$$$ include 'COMMON.GEO'
1700 c$$$ include 'COMMON.MINIM'
1701 c$$$ include 'COMMON.CHAIN'
1703 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1704 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1705 c$$$ + orig_ss_dist(maxres2,maxres2)
1707 c$$$ double precision etot
1708 c$$$ integer iretcode,nfun,i_in,j_in
1711 c$$$ double precision dist
1712 c$$$ external ss_func,fdum
1713 c$$$ double precision ss_func,fdum
1715 c$$$ integer iv(liv),uiparm(2)
1716 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1720 c$$$ call deflt(2,iv,liv,lv,v)
1721 c$$$* 12 means fresh start, dont call deflt
1723 c$$$* max num of fun calls
1724 c$$$ if (maxfun.eq.0) maxfun=500
1726 c$$$* max num of iterations
1727 c$$$ if (maxmin.eq.0) maxmin=1000
1729 c$$$* controls output
1731 c$$$* selects output unit
1734 c$$$* 1 means to print out result
1736 c$$$* 1 means to print out summary stats
1738 c$$$* 1 means to print initial x and d
1740 c$$$* min val for v(radfac) default is 0.1
1742 c$$$* max val for v(radfac) default is 4.0
1745 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1746 c$$$* the sumsl default is 0.1
1748 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1749 c$$$* the sumsl default is 100*machep
1750 c$$$ v(34)=v(34)/100.0D0
1751 c$$$* absolute convergence
1752 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1755 c$$$* relative convergence
1756 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1759 c$$$* controls initial step size
1761 c$$$* large vals of d correspond to small components of step
1768 c$$$ orig_ss_dc(j,i)=dc(j,i)
1771 c$$$ call geom_to_var(nvar,orig_ss_var)
1775 c$$$ orig_ss_dist(j,i)=dist(j,i)
1776 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1777 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1778 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1790 c$$$ if (ialph(i,1).gt.0) then
1793 c$$$ x(k)=dc(j,i+nres)
1800 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1803 c$$$ nfun=iv(6)+iv(30)
1813 c$$$ if (ialph(i,1).gt.0) then
1816 c$$$ dc(j,i+nres)=x(k)
1820 c$$$ call chainbuild_cart
1825 c$$$C-----------------------------------------------------------------------------
1827 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1829 c$$$ include 'DIMENSIONS'
1830 c$$$ include 'COMMON.DERIV'
1831 c$$$ include 'COMMON.IOUNITS'
1832 c$$$ include 'COMMON.VAR'
1833 c$$$ include 'COMMON.CHAIN'
1834 c$$$ include 'COMMON.INTERACT'
1835 c$$$ include 'COMMON.SBRIDGE'
1837 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1838 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1839 c$$$ + orig_ss_dist(maxres2,maxres2)
1842 c$$$ double precision x(maxres6)
1844 c$$$ double precision f
1845 c$$$ integer uiparm(2)
1846 c$$$ real*8 urparm(1)
1847 c$$$ external ufparm
1848 c$$$ double precision ufparm
1851 c$$$ double precision dist
1853 c$$$ integer i,j,k,ss_i,ss_j
1854 c$$$ double precision tempf,var(maxvar)
1869 c$$$ if (ialph(i,1).gt.0) then
1872 c$$$ dc(j,i+nres)=x(k)
1876 c$$$ call chainbuild_cart
1878 c$$$ call geom_to_var(nvar,var)
1880 c$$$c Constraints on all angles
1882 c$$$ tempf=var(i)-orig_ss_var(i)
1883 c$$$ f=f+tempf*tempf
1886 c$$$c Constraints on all distances
1888 c$$$ if (i.gt.1) then
1889 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1890 c$$$ f=f+tempf*tempf
1893 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1894 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1895 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1896 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1897 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1898 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1899 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1900 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1904 c$$$c Constraints for the relevant CYS-CYS
1905 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1906 c$$$ f=f+tempf*tempf
1907 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
1909 c$$$c$$$ if (nf.ne.nfl) then
1910 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1911 c$$$c$$$ + f,dist(5+nres,14+nres)
1919 c$$$C-----------------------------------------------------------------------------
1920 c$$$C-----------------------------------------------------------------------------
1921 subroutine triple_ssbond_ene(resi,resj,resk,eij)
1922 include 'DIMENSIONS'
1923 include 'COMMON.SBRIDGE'
1924 include 'COMMON.CHAIN'
1925 include 'COMMON.DERIV'
1926 include 'COMMON.LOCAL'
1927 include 'COMMON.INTERACT'
1928 include 'COMMON.VAR'
1929 include 'COMMON.IOUNITS'
1930 include 'COMMON.CALC'
1933 C include 'COMMON.MD'
1937 c External functions
1938 double precision h_base
1942 integer resi,resj,resk
1945 double precision eij,eij1,eij2,eij3
1949 c integer itypi,itypj,k,l
1950 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
1951 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
1952 double precision xik,yik,zik,xjk,yjk,zjk
1953 double precision sig0ij,ljd,sig,fac,e1,e2
1954 double precision dcosom1(3),dcosom2(3),ed
1955 double precision pom1,pom2
1956 double precision ljA,ljB,ljXs
1957 double precision d_ljB(1:3)
1958 double precision ssA,ssB,ssC,ssXs
1959 double precision ssxm,ljxm,ssm,ljm
1960 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
1965 C write(iout,*) resi,resj,resk
1967 dxi=dc_norm(1,nres+i)
1968 dyi=dc_norm(2,nres+i)
1969 dzi=dc_norm(3,nres+i)
1970 dsci_inv=vbld_inv(i+nres)
1980 dxj=dc_norm(1,nres+j)
1981 dyj=dc_norm(2,nres+j)
1982 dzj=dc_norm(3,nres+j)
1983 dscj_inv=vbld_inv(j+nres)
1989 dxk=dc_norm(1,nres+k)
1990 dyk=dc_norm(2,nres+k)
1991 dzk=dc_norm(3,nres+k)
1992 dscj_inv=vbld_inv(k+nres)
2002 rrij=(xij*xij+yij*yij+zij*zij)
2003 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2004 rrik=(xik*xik+yik*yik+zik*zik)
2006 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2008 C there are three combination of distances for each trisulfide bonds
2009 C The first case the ith atom is the center
2010 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2011 C distance y is second distance the a,b,c,d are parameters derived for
2012 C this problem d parameter was set as a penalty currenlty set to 1.
2013 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2014 C second case jth atom is center
2015 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2016 C the third case kth atom is the center
2017 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2022 C write(iout,*)i,j,k,eij
2023 C The energy penalty calculated now time for the gradient part
2024 C derivative over rij
2025 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2026 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2031 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2032 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2035 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2036 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2038 C now derivative over rik
2039 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2040 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2045 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2046 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2049 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2050 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2052 C now derivative over rjk
2053 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2054 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2059 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2060 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2063 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2064 gvdwc(l,k)=gvdwc(l,k)+gg(l)