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.SETUP'
94 include 'COMMON.SBRIDGE'
95 include 'COMMON.CHAIN'
96 include 'COMMON.DERIV'
97 include 'COMMON.LOCAL'
98 include 'COMMON.INTERACT'
100 include 'COMMON.IOUNITS'
101 include 'COMMON.CALC'
102 include 'COMMON.NAMES'
110 double precision h_base
121 c integer itypi,itypj,k,l
122 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
123 double precision sig0ij,ljd,sig,fac,e1,e2
124 double precision dcosom1(3),dcosom2(3),ed
125 double precision pom1,pom2
126 double precision ljA,ljB,ljXs
127 double precision d_ljB(1:3)
128 double precision ssA,ssB,ssC,ssXs
129 double precision ssxm,ljxm,ssm,ljm
130 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
131 double precision f1,f2,h1,h2,hd1,hd2
132 double precision omega,delta_inv,deltasq_inv,fac1,fac2
134 double precision xm,d_xm(1:3)
135 double precision sslipi,sslipj,ssgradlipi,ssgradlipj
136 integer ici,icj,itypi,itypj
137 double precision boxshift,sscale,sscagrad
138 c-------END FIRST METHOD
139 c-------SECOND METHOD
140 c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
141 c-------END SECOND METHOD
144 logical checkstop,transgrad
145 common /sschecks/ checkstop,transgrad
147 integer icheck,nicheck,jcheck,njcheck
148 double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
149 c-------END TESTING CODE
156 if (ici.eq.0 .or. icj.eq.0) then
158 write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
159 & "Processor",me," attempt to create",
160 & " a disulfide link between non-cysteine residues ",restyp(i),i,
162 call MPI_Abort(MPI_COMM_WORLD,ierr)
164 write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
165 & "Processor",me," attempt to create",
166 & " a disulfide link between non-cysteine residues ",restyp(i),i,
172 dxi=dc_norm(1,nres+i)
173 dyi=dc_norm(2,nres+i)
174 dzi=dc_norm(3,nres+i)
175 dsci_inv=vbld_inv(i+nres)
179 call to_box(xi,yi,zi)
180 C define scaling factor for lipids
182 C if (positi.le.0) positi=positi+boxzsize
184 C first for peptide groups
185 c for each residue check if it is in lipid or lipid water border area
186 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
191 call to_box(xj,yj,zj)
192 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
193 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
194 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
195 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
196 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
197 xj=boxshift(xj-xi,boxxsize)
198 yj=boxshift(yj-yi,boxysize)
199 zj=boxshift(zj-zi,boxzsize)
200 dxj=dc_norm(1,nres+j)
201 dyj=dc_norm(2,nres+j)
202 dzj=dc_norm(3,nres+j)
203 dscj_inv=vbld_inv(j+nres)
205 chi1=chi(itypi,itypj)
206 chi2=chi(itypj,itypi)
213 alf12=0.5D0*(alf1+alf2)
215 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
216 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
217 sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
218 sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
219 c The following are set in sc_angular
223 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
224 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
225 c om12=dxi*dxj+dyi*dyj+dzi*dzj
227 rij=1.0D0/rij ! Reset this so it makes sense
229 sig0ij=sigma(itypi,itypj)
230 sig=sig0ij*dsqrt(1.0D0/sigsq)
233 ljA=eps1*eps2rt**2*eps3rt**2
236 ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
241 deltat12=om2-om1+2.0d0
246 & +akth*(deltat1*deltat1+deltat2*deltat2)
247 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
248 ssxm=ssXs-0.5D0*ssB/ssA
251 c$$$c Some extra output
252 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
253 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
254 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
255 c$$$ if (ssx0.gt.0.0d0) then
256 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
260 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
261 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
262 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
264 c-------END TESTING CODE
267 c Stop and plot energy and derivative as a function of distance
269 ssm=ssC-0.25D0*ssB*ssB/ssA
270 ljm=-0.25D0*ljB*bb/aa
272 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
280 if (.not.checkstop) then
287 if (checkstop) rij=(ssxm-1.0d0)+
288 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
289 c-------END TESTING CODE
291 if (rij.gt.ljxm) then
294 fac=(1.0D0/ljd)**expon
297 eij=eps1*eps2rt*eps3rt*(e1+e2)
300 eij=eij*eps2rt*eps3rt*sss
303 e1=e1*eps1*eps2rt**2*eps3rt**2
304 ed=-expon*(e1+eij)/ljd
306 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
307 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
308 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
309 eom12=eij*eps1_om12+eps2der*eps2rt_om12
310 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
311 else if (rij.lt.ssxm) then
314 eij=ssA*ssd*ssd+ssB*ssd+ssC
316 ed=2*akcm*ssd+akct*deltat12
317 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
319 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
320 eom1=-2*akth*deltat1-pom1-om2*pom2
321 eom2= 2*akth*deltat2+pom1-om1*pom2
324 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
326 d_ssxm(1)=0.5D0*akct/ssA
330 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
331 d_ljxm(2)=d_ljxm(1)*sigsq_om2
332 d_ljxm(3)=d_ljxm(1)*sigsq_om12
333 d_ljxm(1)=d_ljxm(1)*sigsq_om1
335 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
338 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
342 ssm=ssC-0.25D0*ssB*ssB/ssA
343 d_ssm(1)=0.5D0*akct*ssB/ssA
344 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
345 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
347 f1=(rij-xm)/(ssxm-xm)
348 f2=(rij-ssxm)/(xm-ssxm)
352 delta_inv=1.0d0/(xm-ssxm)
353 deltasq_inv=delta_inv*delta_inv
355 fac1=deltasq_inv*fac*(xm-rij)
356 fac2=deltasq_inv*fac*(rij-ssxm)
357 ed=delta_inv*(Ht*hd2-ssm*hd1)
359 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
360 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
361 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
362 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
365 ljm=-0.25D0*ljB*bb/aa
366 d_ljm(1)=-0.5D0*bb/aa*ljB
367 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
368 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
370 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
371 f1=(rij-ljxm)/(xm-ljxm)
372 f2=(rij-xm)/(ljxm-xm)
376 delta_inv=1.0d0/(ljxm-xm)
377 deltasq_inv=delta_inv*delta_inv
379 fac1=deltasq_inv*fac*(ljxm-rij)
380 fac2=deltasq_inv*fac*(rij-xm)
381 ed=delta_inv*(ljm*hd2-Ht*hd1)
383 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
384 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
385 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
386 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
388 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
390 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
396 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
397 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
398 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
400 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
401 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
402 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
403 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
406 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
408 c$$$ d_ljm(k)=ljm*d_ljB(k)
412 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
413 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
414 c$$$ d_ss(2)=akct*ssd
415 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
416 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
419 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
420 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
421 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
423 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
424 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
426 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
428 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
429 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
430 c$$$ h1=h_base(f1,hd1)
431 c$$$ h2=h_base(f2,hd2)
432 c$$$ eij=ss*h1+ljf*h2
433 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
434 c$$$ deltasq_inv=delta_inv*delta_inv
435 c$$$ fac=ljf*hd2-ss*hd1
436 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
437 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
438 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
439 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
440 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
441 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
442 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
444 c$$$ havebond=.false.
445 c$$$ if (ed.gt.0.0d0) havebond=.true.
446 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
453 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
454 c write(iout,'(a15,f12.2,f8.1,2i5)')
455 c & "SSBOND_E_FORM",totT,t_bath,i,j
459 dyn_ssbond_ij(ici,icj)=eij
460 else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300)
462 dyn_ssbond_ij(ici,icj)=1.0d300
465 c write(iout,'(a15,f12.2,f8.1,2i5)')
466 c & "SSBOND_E_BREAK",totT,t_bath,i,j
473 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
474 & "CHECKSTOP",rij,eij,ed
479 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
486 c-------END TESTING CODE
487 gg_lipi(3)=ssgradlipi*eij
488 gg_lipj(3)=ssgradlipj*eij
491 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
492 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
495 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
498 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
499 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
500 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
501 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
502 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
503 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
507 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
512 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
513 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
519 C-----------------------------------------------------------------------------
521 double precision function h_base(x,deriv)
522 c A smooth function going 0->1 in range [0,1]
523 c It should NOT be called outside range [0,1], it will not work there.
530 double precision deriv
536 c Two parabolas put together. First derivative zero at extrema
537 c$$$ if (x.lt.0.5D0) then
538 c$$$ h_base=2.0D0*x*x
542 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
543 c$$$ deriv=4.0D0*deriv
546 c Third degree polynomial. First derivative zero at extrema
547 h_base=x*x*(3.0d0-2.0d0*x)
548 deriv=6.0d0*x*(1.0d0-x)
550 c Fifth degree polynomial. First and second derivatives zero at extrema
552 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
554 c$$$ deriv=deriv*deriv
555 c$$$ deriv=30.0d0*xsq*deriv
560 c----------------------------------------------------------------------------
562 subroutine dyn_set_nss
563 c Adjust nss and other relevant variables based on dyn_ssbond_ij
571 include 'COMMON.SBRIDGE'
572 include 'COMMON.CHAIN'
573 include 'COMMON.IOUNITS'
574 include 'COMMON.SETUP'
582 double precision emin
584 integer diff,allflag(maxss),allnss,
585 & allihpb(maxss),alljhpb(maxss),
586 & newnss,newihpb(maxss),newjhpb(maxss)
588 integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
589 integer g_newihpb(maxss),g_newjhpb(maxss),g_newnss
594 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
603 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
607 if (allflag(i).eq.0 .and.
608 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
609 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
613 if (emin.lt.1.0d300) then
616 if (allflag(i).eq.0 .and.
617 & (allihpb(i).eq.allihpb(imin) .or.
618 & alljhpb(i).eq.allihpb(imin) .or.
619 & allihpb(i).eq.alljhpb(imin) .or.
620 & alljhpb(i).eq.alljhpb(imin))) then
627 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
631 if (allflag(i).eq.1) then
633 newihpb(newnss)=allihpb(i)
634 newjhpb(newnss)=alljhpb(i)
639 if (nfgtasks.gt.1)then
641 call MPI_Reduce(newnss,g_newnss,1,
642 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
643 call MPI_Gather(newnss,1,MPI_INTEGER,
644 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
647 displ(i)=i_newnss(i-1)+displ(i-1)
649 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
650 & g_newihpb,i_newnss,displ,MPI_INTEGER,
652 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
653 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
655 if(fg_rank.eq.0) then
656 c print *,'g_newnss',g_newnss
657 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
658 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
661 newihpb(i)=g_newihpb(i)
662 newjhpb(i)=g_newjhpb(i)
670 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
675 if (idssb(i).eq.newihpb(j) .and.
676 & jdssb(i).eq.newjhpb(j)) found=.true.
680 c if (.not.found.and.fg_rank.eq.0)
681 c & write(iout,'(a15,f12.2,f8.1,2i5)')
682 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
690 if (newihpb(i).eq.idssb(j) .and.
691 & newjhpb(i).eq.jdssb(j)) found=.true.
695 c if (.not.found.and.fg_rank.eq.0)
696 c & write(iout,'(a15,f12.2,f8.1,2i5)')
697 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
712 C-----------------------------------------------------------------------------
713 C-----------------------------------------------------------------------------
714 C-----------------------------------------------------------------------------
716 c$$$c-----------------------------------------------------------------------------
718 c$$$ subroutine ss_relax(i_in,j_in)
722 c$$$ include 'DIMENSIONS'
723 c$$$ include 'COMMON.VAR'
724 c$$$ include 'COMMON.CHAIN'
725 c$$$ include 'COMMON.IOUNITS'
726 c$$$ include 'COMMON.INTERACT'
728 c$$$c Input arguments
729 c$$$ integer i_in,j_in
731 c$$$c Local variables
732 c$$$ integer i,iretcode,nfun_sc
734 c$$$ double precision var(maxvar),e_sc,etot
741 c$$$ mask_side(i_in)=1
742 c$$$ mask_side(j_in)=1
744 c$$$c Minimize the two selected side-chains
745 c$$$ call overlap_sc(scfail) ! Better not fail!
746 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
753 c$$$c-------------------------------------------------------------
755 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
756 c$$$c Minimize side-chains only, starting from geom but without modifying
758 c$$$c If mask_r is already set, only the selected side-chains are minimized,
759 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
763 c$$$ include 'DIMENSIONS'
764 c$$$ include 'COMMON.IOUNITS'
765 c$$$ include 'COMMON.VAR'
766 c$$$ include 'COMMON.CHAIN'
767 c$$$ include 'COMMON.GEO'
768 c$$$ include 'COMMON.MINIM'
770 c$$$ common /srutu/ icall
772 c$$$c Output arguments
773 c$$$ double precision etot_sc
774 c$$$ integer iretcode,nfun
776 c$$$c External functions/subroutines
777 c$$$ external func_sc,grad_sc,fdum
779 c$$$c Local variables
781 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
783 c$$$ double precision rdum(1)
784 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
786 c$$$ integer i,nvar_restr
789 c$$$cmc start_minim=.true.
790 c$$$ call deflt(2,iv,liv,lv,v)
791 c$$$* 12 means fresh start, dont call deflt
793 c$$$* max num of fun calls
794 c$$$ if (maxfun.eq.0) maxfun=500
796 c$$$* max num of iterations
797 c$$$ if (maxmin.eq.0) maxmin=1000
799 c$$$* controls output
801 c$$$* selects output unit
803 c$$$c iv(21)=iout ! DEBUG
804 c$$$c iv(21)=8 ! DEBUG
805 c$$$* 1 means to print out result
807 c$$$c iv(22)=1 ! DEBUG
808 c$$$* 1 means to print out summary stats
810 c$$$c iv(23)=1 ! DEBUG
811 c$$$* 1 means to print initial x and d
813 c$$$c iv(24)=1 ! DEBUG
814 c$$$* min val for v(radfac) default is 0.1
816 c$$$* max val for v(radfac) default is 4.0
819 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
820 c$$$* the sumsl default is 0.1
822 c$$$* false conv if (act fnctn decrease) .lt. v(34)
823 c$$$* the sumsl default is 100*machep
824 c$$$ v(34)=v(34)/100.0D0
825 c$$$* absolute convergence
826 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
828 c$$$* relative convergence
829 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
831 c$$$* controls initial step size
833 c$$$* large vals of d correspond to small components of step
837 c$$$ do i=nphi+1,nvar
841 c$$$ call geom_to_var(nvar,x)
842 c$$$ IF (mask_r) THEN
843 c$$$ do i=1,nres ! Just in case...
847 c$$$ call x2xx(x,xx,nvar_restr)
848 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
849 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
852 c$$$c When minimizing ALL side-chains, etotal_sc is a little
853 c$$$c faster if we don't set mask_r
859 c$$$ call x2xx(x,xx,nvar_restr)
860 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
861 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
864 c$$$ call var_to_geom(nvar,x)
865 c$$$ call chainbuild_sc
872 c$$$C--------------------------------------------------------------------------
874 c$$$ subroutine chainbuild_sc
876 c$$$ include 'DIMENSIONS'
877 c$$$ include 'COMMON.VAR'
878 c$$$ include 'COMMON.INTERACT'
880 c$$$c Local variables
885 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
886 c$$$ call locate_side_chain(i)
893 c$$$C--------------------------------------------------------------------------
895 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
899 c$$$ include 'DIMENSIONS'
900 c$$$ include 'COMMON.DERIV'
901 c$$$ include 'COMMON.VAR'
902 c$$$ include 'COMMON.MINIM'
903 c$$$ include 'COMMON.IOUNITS'
905 c$$$c Input arguments
907 c$$$ double precision x(maxvar)
908 c$$$ double precision ufparm
911 c$$$c Input/Output arguments
913 c$$$ integer uiparm(1)
914 c$$$ double precision urparm(1)
916 c$$$c Output arguments
917 c$$$ double precision f
919 c$$$c Local variables
920 c$$$ double precision energia(0:n_ene)
922 c$$$c Variables used to intercept NaNs
923 c$$$ double precision x_sum
932 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
935 c$$$ x_sum=x_sum+x(i_NAN)
937 c$$$c Calculate the energy only if the coordinates are ok
938 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
939 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
945 c$$$ call var_to_geom_restr(n,x)
947 c$$$ call chainbuild_sc
948 c$$$ call etotal_sc(energia(0))
950 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
959 c$$$c-------------------------------------------------------
961 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
965 c$$$ include 'DIMENSIONS'
966 c$$$ include 'COMMON.CHAIN'
967 c$$$ include 'COMMON.DERIV'
968 c$$$ include 'COMMON.VAR'
969 c$$$ include 'COMMON.INTERACT'
970 c$$$ include 'COMMON.MINIM'
972 c$$$c Input arguments
974 c$$$ double precision x(maxvar)
975 c$$$ double precision ufparm
978 c$$$c Input/Output arguments
980 c$$$ integer uiparm(1)
981 c$$$ double precision urparm(1)
983 c$$$c Output arguments
984 c$$$ double precision g(maxvar)
986 c$$$c Local variables
987 c$$$ double precision f,gphii,gthetai,galphai,gomegai
988 c$$$ integer ig,ind,i,j,k,igall,ij
992 c$$$ if (nf-nfl+1) 20,30,40
993 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
994 c$$$c write (iout,*) 'grad 20'
995 c$$$ if (nf.eq.0) return
997 c$$$ 30 call var_to_geom_restr(n,x)
998 c$$$ call chainbuild_sc
1000 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1002 c$$$ 40 call cartder
1004 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1010 c$$$ IF (mask_phi(i+2).eq.1) THEN
1012 c$$$ do j=i+1,nres-1
1015 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1016 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1022 c$$$ ind=ind+nres-1-i
1029 c$$$ IF (mask_theta(i+2).eq.1) THEN
1032 c$$$ do j=i+1,nres-1
1035 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1036 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1041 c$$$ ind=ind+nres-1-i
1046 c$$$ if (itype(i).ne.10) then
1047 c$$$ IF (mask_side(i).eq.1) THEN
1051 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1060 c$$$ if (itype(i).ne.10) then
1061 c$$$ IF (mask_side(i).eq.1) THEN
1065 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1073 c$$$C Add the components corresponding to local energy terms.
1080 c$$$ if (mask_phi(i).eq.1) then
1082 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1088 c$$$ if (mask_theta(i).eq.1) then
1090 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1096 c$$$ if (itype(i).ne.10) then
1098 c$$$ if (mask_side(i).eq.1) then
1100 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1107 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1113 c$$$C-----------------------------------------------------------------------------
1115 c$$$ subroutine etotal_sc(energy_sc)
1119 c$$$ include 'DIMENSIONS'
1120 c$$$ include 'COMMON.VAR'
1121 c$$$ include 'COMMON.INTERACT'
1122 c$$$ include 'COMMON.DERIV'
1123 c$$$ include 'COMMON.FFIELD'
1125 c$$$c Output arguments
1126 c$$$ double precision energy_sc(0:n_ene)
1128 c$$$c Local variables
1129 c$$$ double precision evdw,escloc
1134 c$$$ energy_sc(i)=0.0D0
1137 c$$$ if (mask_r) then
1138 c$$$ call egb_sc(evdw)
1139 c$$$ call esc_sc(escloc)
1142 c$$$ call esc(escloc)
1145 c$$$ if (evdw.eq.1.0D20) then
1146 c$$$ energy_sc(0)=evdw
1148 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1150 c$$$ energy_sc(1)=evdw
1151 c$$$ energy_sc(12)=escloc
1154 c$$$C Sum up the components of the Cartesian gradient.
1158 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1165 c$$$C-----------------------------------------------------------------------------
1167 c$$$ subroutine egb_sc(evdw)
1169 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1170 c$$$C assuming the Gay-Berne potential of interaction.
1172 c$$$ implicit real*8 (a-h,o-z)
1173 c$$$ include 'DIMENSIONS'
1174 c$$$ include 'COMMON.GEO'
1175 c$$$ include 'COMMON.VAR'
1176 c$$$ include 'COMMON.LOCAL'
1177 c$$$ include 'COMMON.CHAIN'
1178 c$$$ include 'COMMON.DERIV'
1179 c$$$ include 'COMMON.NAMES'
1180 c$$$ include 'COMMON.INTERACT'
1181 c$$$ include 'COMMON.IOUNITS'
1182 c$$$ include 'COMMON.CALC'
1183 c$$$ include 'COMMON.CONTROL'
1186 c$$$ energy_dec=.false.
1187 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1190 c$$$c if (icall.eq.0) lprn=.false.
1192 c$$$ do i=iatsc_s,iatsc_e
1194 c$$$ itypi1=itype(i+1)
1198 c$$$ dxi=dc_norm(1,nres+i)
1199 c$$$ dyi=dc_norm(2,nres+i)
1200 c$$$ dzi=dc_norm(3,nres+i)
1201 c$$$c dsci_inv=dsc_inv(itypi)
1202 c$$$ dsci_inv=vbld_inv(i+nres)
1203 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1204 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1206 c$$$C Calculate SC interaction energy.
1208 c$$$ do iint=1,nint_gr(i)
1209 c$$$ do j=istart(i,iint),iend(i,iint)
1210 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1213 c$$$c dscj_inv=dsc_inv(itypj)
1214 c$$$ dscj_inv=vbld_inv(j+nres)
1215 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1216 c$$$c & 1.0d0/vbld(j+nres)
1217 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1218 c$$$ sig0ij=sigma(itypi,itypj)
1219 c$$$ chi1=chi(itypi,itypj)
1220 c$$$ chi2=chi(itypj,itypi)
1221 c$$$ chi12=chi1*chi2
1222 c$$$ chip1=chip(itypi)
1223 c$$$ chip2=chip(itypj)
1224 c$$$ chip12=chip1*chip2
1225 c$$$ alf1=alp(itypi)
1226 c$$$ alf2=alp(itypj)
1227 c$$$ alf12=0.5D0*(alf1+alf2)
1228 c$$$C For diagnostics only!!!
1238 c$$$ xj=c(1,nres+j)-xi
1239 c$$$ yj=c(2,nres+j)-yi
1240 c$$$ zj=c(3,nres+j)-zi
1241 c$$$ dxj=dc_norm(1,nres+j)
1242 c$$$ dyj=dc_norm(2,nres+j)
1243 c$$$ dzj=dc_norm(3,nres+j)
1244 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1245 c$$$c write (iout,*) "j",j," dc_norm",
1246 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1247 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1248 c$$$ rij=dsqrt(rrij)
1249 c$$$C Calculate angle-dependent terms of energy and contributions to their
1251 c$$$ call sc_angular
1252 c$$$ sigsq=1.0D0/sigsq
1253 c$$$ sig=sig0ij*dsqrt(sigsq)
1254 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1255 c$$$c for diagnostics; uncomment
1256 c$$$c rij_shift=1.2*sig0ij
1257 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1258 c$$$ if (rij_shift.le.0.0D0) then
1260 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1261 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1262 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1265 c$$$ sigder=-sig*sigsq
1266 c$$$c---------------------------------------------------------------
1267 c$$$ rij_shift=1.0D0/rij_shift
1268 c$$$ fac=rij_shift**expon
1269 c$$$ e1=fac*fac*aa(itypi,itypj)
1270 c$$$ e2=fac*bb(itypi,itypj)
1271 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1272 c$$$ eps2der=evdwij*eps3rt
1273 c$$$ eps3der=evdwij*eps2rt
1274 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1275 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1276 c$$$ evdwij=evdwij*eps2rt*eps3rt
1277 c$$$ evdw=evdw+evdwij
1279 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1280 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1281 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1282 c$$$ & restyp(itypi),i,restyp(itypj),j,
1283 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1284 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1285 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1289 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1290 c$$$ & 'evdw',i,j,evdwij
1292 c$$$C Calculate gradient components.
1293 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1294 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1295 c$$$ sigder=fac*sigder
1298 c$$$C Calculate the radial part of the gradient
1302 c$$$C Calculate angular part of the gradient.
1308 c$$$ energy_dec=.false.
1312 c$$$c-----------------------------------------------------------------------------
1314 c$$$ subroutine esc_sc(escloc)
1315 c$$$C Calculate the local energy of a side chain and its derivatives in the
1316 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1317 c$$$C ALPHA and OMEGA.
1318 c$$$ implicit real*8 (a-h,o-z)
1319 c$$$ include 'DIMENSIONS'
1320 c$$$ include 'COMMON.GEO'
1321 c$$$ include 'COMMON.LOCAL'
1322 c$$$ include 'COMMON.VAR'
1323 c$$$ include 'COMMON.INTERACT'
1324 c$$$ include 'COMMON.DERIV'
1325 c$$$ include 'COMMON.CHAIN'
1326 c$$$ include 'COMMON.IOUNITS'
1327 c$$$ include 'COMMON.NAMES'
1328 c$$$ include 'COMMON.FFIELD'
1329 c$$$ include 'COMMON.CONTROL'
1330 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1331 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1332 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1333 c$$$ delta=0.02d0*pi
1335 c$$$c write (iout,'(a)') 'ESC'
1336 c$$$ do i=loc_start,loc_end
1337 c$$$ IF (mask_side(i).eq.1) THEN
1339 c$$$ if (it.eq.10) goto 1
1340 c$$$ nlobit=nlob(it)
1341 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1342 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1343 c$$$ theti=theta(i+1)-pipol
1344 c$$$ x(1)=dtan(theti)
1348 c$$$ if (x(2).gt.pi-delta) then
1350 c$$$ xtemp(2)=pi-delta
1352 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1354 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1355 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1356 c$$$ & escloci,dersc(2))
1357 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1358 c$$$ & ddersc0(1),dersc(1))
1359 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1360 c$$$ & ddersc0(3),dersc(3))
1361 c$$$ xtemp(2)=pi-delta
1362 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1364 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1365 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1366 c$$$ & dersc0(2),esclocbi,dersc02)
1367 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1368 c$$$ & dersc12,dersc01)
1369 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1370 c$$$ dersc0(1)=dersc01
1371 c$$$ dersc0(2)=dersc02
1372 c$$$ dersc0(3)=0.0d0
1374 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1376 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1377 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1378 c$$$c & esclocbi,ss,ssd
1379 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1380 c$$$c escloci=esclocbi
1381 c$$$c write (iout,*) escloci
1382 c$$$ else if (x(2).lt.delta) then
1386 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1388 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1389 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1390 c$$$ & escloci,dersc(2))
1391 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1392 c$$$ & ddersc0(1),dersc(1))
1393 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1394 c$$$ & ddersc0(3),dersc(3))
1396 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1398 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1399 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1400 c$$$ & dersc0(2),esclocbi,dersc02)
1401 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1402 c$$$ & dersc12,dersc01)
1403 c$$$ dersc0(1)=dersc01
1404 c$$$ dersc0(2)=dersc02
1405 c$$$ dersc0(3)=0.0d0
1406 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1408 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1410 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1411 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1412 c$$$c & esclocbi,ss,ssd
1413 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1414 c$$$c write (iout,*) escloci
1416 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1419 c$$$ escloc=escloc+escloci
1420 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1421 c$$$ & 'escloc',i,escloci
1422 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1424 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1425 c$$$ & wscloc*dersc(1)
1426 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1427 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1434 c$$$C-----------------------------------------------------------------------------
1436 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1438 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1439 c$$$C assuming the Gay-Berne potential of interaction.
1441 c$$$ implicit real*8 (a-h,o-z)
1442 c$$$ include 'DIMENSIONS'
1443 c$$$ include 'COMMON.GEO'
1444 c$$$ include 'COMMON.VAR'
1445 c$$$ include 'COMMON.LOCAL'
1446 c$$$ include 'COMMON.CHAIN'
1447 c$$$ include 'COMMON.DERIV'
1448 c$$$ include 'COMMON.NAMES'
1449 c$$$ include 'COMMON.INTERACT'
1450 c$$$ include 'COMMON.IOUNITS'
1451 c$$$ include 'COMMON.CALC'
1452 c$$$ include 'COMMON.CONTROL'
1455 c$$$ energy_dec=.false.
1456 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1460 c$$$c$$$ do i=iatsc_s,iatsc_e
1463 c$$$ itypi1=itype(i+1)
1467 c$$$ dxi=dc_norm(1,nres+i)
1468 c$$$ dyi=dc_norm(2,nres+i)
1469 c$$$ dzi=dc_norm(3,nres+i)
1470 c$$$c dsci_inv=dsc_inv(itypi)
1471 c$$$ dsci_inv=vbld_inv(i+nres)
1472 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1473 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1475 c$$$C Calculate SC interaction energy.
1477 c$$$c$$$ do iint=1,nint_gr(i)
1478 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1482 c$$$c dscj_inv=dsc_inv(itypj)
1483 c$$$ dscj_inv=vbld_inv(j+nres)
1484 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1485 c$$$c & 1.0d0/vbld(j+nres)
1486 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1487 c$$$ sig0ij=sigma(itypi,itypj)
1488 c$$$ chi1=chi(itypi,itypj)
1489 c$$$ chi2=chi(itypj,itypi)
1490 c$$$ chi12=chi1*chi2
1491 c$$$ chip1=chip(itypi)
1492 c$$$ chip2=chip(itypj)
1493 c$$$ chip12=chip1*chip2
1494 c$$$ alf1=alp(itypi)
1495 c$$$ alf2=alp(itypj)
1496 c$$$ alf12=0.5D0*(alf1+alf2)
1497 c$$$C For diagnostics only!!!
1507 c$$$ xj=c(1,nres+j)-xi
1508 c$$$ yj=c(2,nres+j)-yi
1509 c$$$ zj=c(3,nres+j)-zi
1510 c$$$ dxj=dc_norm(1,nres+j)
1511 c$$$ dyj=dc_norm(2,nres+j)
1512 c$$$ dzj=dc_norm(3,nres+j)
1513 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1514 c$$$c write (iout,*) "j",j," dc_norm",
1515 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1516 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1517 c$$$ rij=dsqrt(rrij)
1518 c$$$C Calculate angle-dependent terms of energy and contributions to their
1520 c$$$ call sc_angular
1521 c$$$ sigsq=1.0D0/sigsq
1522 c$$$ sig=sig0ij*dsqrt(sigsq)
1523 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1524 c$$$c for diagnostics; uncomment
1525 c$$$c rij_shift=1.2*sig0ij
1526 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1527 c$$$ if (rij_shift.le.0.0D0) then
1529 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1530 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1531 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1534 c$$$ sigder=-sig*sigsq
1535 c$$$c---------------------------------------------------------------
1536 c$$$ rij_shift=1.0D0/rij_shift
1537 c$$$ fac=rij_shift**expon
1538 c$$$ e1=fac*fac*aa(itypi,itypj)
1539 c$$$ e2=fac*bb(itypi,itypj)
1540 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1541 c$$$ eps2der=evdwij*eps3rt
1542 c$$$ eps3der=evdwij*eps2rt
1543 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1544 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1545 c$$$ evdwij=evdwij*eps2rt*eps3rt
1546 c$$$ evdw=evdw+evdwij
1548 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1549 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1550 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1551 c$$$ & restyp(itypi),i,restyp(itypj),j,
1552 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1553 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1554 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1558 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1559 c$$$ & 'evdw',i,j,evdwij
1561 c$$$C Calculate gradient components.
1562 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1563 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1564 c$$$ sigder=fac*sigder
1567 c$$$C Calculate the radial part of the gradient
1571 c$$$C Calculate angular part of the gradient.
1574 c$$$c$$$ enddo ! iint
1576 c$$$ energy_dec=.false.
1580 c$$$C-----------------------------------------------------------------------------
1582 c$$$ subroutine perturb_side_chain(i,angle)
1586 c$$$ include 'DIMENSIONS'
1587 c$$$ include 'COMMON.CHAIN'
1588 c$$$ include 'COMMON.GEO'
1589 c$$$ include 'COMMON.VAR'
1590 c$$$ include 'COMMON.LOCAL'
1591 c$$$ include 'COMMON.IOUNITS'
1593 c$$$c External functions
1594 c$$$ external ran_number
1595 c$$$ double precision ran_number
1597 c$$$c Input arguments
1599 c$$$ double precision angle ! In degrees
1601 c$$$c Local variables
1603 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1607 c$$$ rad_ang=angle*deg2rad
1610 c$$$ do while (length.lt.0.01)
1611 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1612 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1613 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1614 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1615 c$$$ + rand_v(3)*rand_v(3)
1616 c$$$ length=sqrt(length)
1617 c$$$ rand_v(1)=rand_v(1)/length
1618 c$$$ rand_v(2)=rand_v(2)/length
1619 c$$$ rand_v(3)=rand_v(3)/length
1620 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1621 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1622 c$$$ length=1.0D0-cost*cost
1623 c$$$ if (length.lt.0.0D0) length=0.0D0
1624 c$$$ length=sqrt(length)
1625 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1626 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1627 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1629 c$$$ rand_v(1)=rand_v(1)/length
1630 c$$$ rand_v(2)=rand_v(2)/length
1631 c$$$ rand_v(3)=rand_v(3)/length
1633 c$$$ cost=dcos(rad_ang)
1634 c$$$ sint=dsin(rad_ang)
1635 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1636 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1637 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1638 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1639 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1640 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1641 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1642 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1643 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1645 c$$$ call chainbuild_cart
1650 c$$$c----------------------------------------------------------------------------
1652 c$$$ subroutine ss_relax3(i_in,j_in)
1656 c$$$ include 'DIMENSIONS'
1657 c$$$ include 'COMMON.VAR'
1658 c$$$ include 'COMMON.CHAIN'
1659 c$$$ include 'COMMON.IOUNITS'
1660 c$$$ include 'COMMON.INTERACT'
1662 c$$$c External functions
1663 c$$$ external ran_number
1664 c$$$ double precision ran_number
1666 c$$$c Input arguments
1667 c$$$ integer i_in,j_in
1669 c$$$c Local variables
1670 c$$$ double precision energy_sc(0:n_ene),etot
1671 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1672 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1673 c$$$ integer n,i_pert,i
1674 c$$$ logical notdone
1683 c$$$ mask_side(i_in)=1
1684 c$$$ mask_side(j_in)=1
1686 c$$$ call etotal_sc(energy_sc)
1687 c$$$ etot=energy_sc(0)
1688 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1689 c$$$c + energy_sc(1),energy_sc(12)
1693 c$$$ do while (notdone)
1694 c$$$ if (mod(n,2).eq.0) then
1702 c$$$ org_dc(i)=dc(i,i_pert+nres)
1703 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1704 c$$$ org_c(i)=c(i,i_pert+nres)
1706 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1707 c$$$ call perturb_side_chain(i_pert,ang_pert)
1708 c$$$ call etotal_sc(energy_sc)
1709 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1710 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1711 c$$$ if (rand_fact.lt.exp_fact) then
1712 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1713 c$$$c + energy_sc(1),energy_sc(12)
1714 c$$$ etot=energy_sc(0)
1716 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1717 c$$$c + energy_sc(1),energy_sc(12)
1719 c$$$ dc(i,i_pert+nres)=org_dc(i)
1720 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1721 c$$$ c(i,i_pert+nres)=org_c(i)
1725 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1733 c$$$c----------------------------------------------------------------------------
1735 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1737 c$$$ include 'DIMENSIONS'
1739 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1740 c$$$*********************************************************************
1741 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1742 c$$$* the calling subprogram. *
1743 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1744 c$$$* calculated in the usual pythagorean way. *
1745 c$$$* absolute convergence occurs when the function is within v(31) of *
1746 c$$$* zero. unless you know the minimum value in advance, abs convg *
1747 c$$$* is probably not useful. *
1748 c$$$* relative convergence is when the model predicts that the function *
1749 c$$$* will decrease by less than v(32)*abs(fun). *
1750 c$$$*********************************************************************
1751 c$$$ include 'COMMON.IOUNITS'
1752 c$$$ include 'COMMON.VAR'
1753 c$$$ include 'COMMON.GEO'
1754 c$$$ include 'COMMON.MINIM'
1755 c$$$ include 'COMMON.CHAIN'
1757 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1758 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1759 c$$$ + orig_ss_dist(maxres2,maxres2)
1761 c$$$ double precision etot
1762 c$$$ integer iretcode,nfun,i_in,j_in
1765 c$$$ double precision dist
1766 c$$$ external ss_func,fdum
1767 c$$$ double precision ss_func,fdum
1769 c$$$ integer iv(liv),uiparm(2)
1770 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1774 c$$$ call deflt(2,iv,liv,lv,v)
1775 c$$$* 12 means fresh start, dont call deflt
1777 c$$$* max num of fun calls
1778 c$$$ if (maxfun.eq.0) maxfun=500
1780 c$$$* max num of iterations
1781 c$$$ if (maxmin.eq.0) maxmin=1000
1783 c$$$* controls output
1785 c$$$* selects output unit
1788 c$$$* 1 means to print out result
1790 c$$$* 1 means to print out summary stats
1792 c$$$* 1 means to print initial x and d
1794 c$$$* min val for v(radfac) default is 0.1
1796 c$$$* max val for v(radfac) default is 4.0
1799 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1800 c$$$* the sumsl default is 0.1
1802 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1803 c$$$* the sumsl default is 100*machep
1804 c$$$ v(34)=v(34)/100.0D0
1805 c$$$* absolute convergence
1806 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1809 c$$$* relative convergence
1810 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1813 c$$$* controls initial step size
1815 c$$$* large vals of d correspond to small components of step
1822 c$$$ orig_ss_dc(j,i)=dc(j,i)
1825 c$$$ call geom_to_var(nvar,orig_ss_var)
1829 c$$$ orig_ss_dist(j,i)=dist(j,i)
1830 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1831 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1832 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1844 c$$$ if (ialph(i,1).gt.0) then
1847 c$$$ x(k)=dc(j,i+nres)
1854 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1857 c$$$ nfun=iv(6)+iv(30)
1867 c$$$ if (ialph(i,1).gt.0) then
1870 c$$$ dc(j,i+nres)=x(k)
1874 c$$$ call chainbuild_cart
1879 c$$$C-----------------------------------------------------------------------------
1881 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1883 c$$$ include 'DIMENSIONS'
1884 c$$$ include 'COMMON.DERIV'
1885 c$$$ include 'COMMON.IOUNITS'
1886 c$$$ include 'COMMON.VAR'
1887 c$$$ include 'COMMON.CHAIN'
1888 c$$$ include 'COMMON.INTERACT'
1889 c$$$ include 'COMMON.SBRIDGE'
1891 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1892 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1893 c$$$ + orig_ss_dist(maxres2,maxres2)
1896 c$$$ double precision x(maxres6)
1898 c$$$ double precision f
1899 c$$$ integer uiparm(2)
1900 c$$$ real*8 urparm(1)
1901 c$$$ external ufparm
1902 c$$$ double precision ufparm
1905 c$$$ double precision dist
1907 c$$$ integer i,j,k,ss_i,ss_j
1908 c$$$ double precision tempf,var(maxvar)
1923 c$$$ if (ialph(i,1).gt.0) then
1926 c$$$ dc(j,i+nres)=x(k)
1930 c$$$ call chainbuild_cart
1932 c$$$ call geom_to_var(nvar,var)
1934 c$$$c Constraints on all angles
1936 c$$$ tempf=var(i)-orig_ss_var(i)
1937 c$$$ f=f+tempf*tempf
1940 c$$$c Constraints on all distances
1942 c$$$ if (i.gt.1) then
1943 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1944 c$$$ f=f+tempf*tempf
1947 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1948 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1949 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1950 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1951 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1952 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1953 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1954 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1958 c$$$c Constraints for the relevant CYS-CYS
1959 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1960 c$$$ f=f+tempf*tempf
1961 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
1963 c$$$c$$$ if (nf.ne.nfl) then
1964 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1965 c$$$c$$$ + f,dist(5+nres,14+nres)
1973 c$$$C-----------------------------------------------------------------------------
1974 c$$$C-----------------------------------------------------------------------------
1975 subroutine triple_ssbond_ene(resi,resj,resk,eij)
1976 include 'DIMENSIONS'
1977 include 'COMMON.SBRIDGE'
1978 include 'COMMON.CHAIN'
1979 include 'COMMON.DERIV'
1980 include 'COMMON.LOCAL'
1981 include 'COMMON.INTERACT'
1982 include 'COMMON.VAR'
1983 include 'COMMON.IOUNITS'
1984 include 'COMMON.CALC'
1991 c External functions
1992 double precision h_base
1996 integer resi,resj,resk
1999 double precision eij,eij1,eij2,eij3
2003 c integer itypi,itypj,k,l
2004 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2005 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2006 double precision xik,yik,zik,xjk,yjk,zjk
2007 double precision sig0ij,ljd,sig,fac,e1,e2
2008 double precision dcosom1(3),dcosom2(3),ed
2009 double precision pom1,pom2
2010 double precision ljA,ljB,ljXs
2011 double precision d_ljB(1:3)
2012 double precision ssA,ssB,ssC,ssXs
2013 double precision ssxm,ljxm,ssm,ljm
2014 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2015 if (dtriss.eq.0) return
2019 C write(iout,*) resi,resj,resk
2021 dxi=dc_norm(1,nres+i)
2022 dyi=dc_norm(2,nres+i)
2023 dzi=dc_norm(3,nres+i)
2024 dsci_inv=vbld_inv(i+nres)
2034 dxj=dc_norm(1,nres+j)
2035 dyj=dc_norm(2,nres+j)
2036 dzj=dc_norm(3,nres+j)
2037 dscj_inv=vbld_inv(j+nres)
2043 dxk=dc_norm(1,nres+k)
2044 dyk=dc_norm(2,nres+k)
2045 dzk=dc_norm(3,nres+k)
2046 dscj_inv=vbld_inv(k+nres)
2056 rrij=(xij*xij+yij*yij+zij*zij)
2057 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2058 rrik=(xik*xik+yik*yik+zik*zik)
2060 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2062 C there are three combination of distances for each trisulfide bonds
2063 C The first case the ith atom is the center
2064 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2065 C distance y is second distance the a,b,c,d are parameters derived for
2066 C this problem d parameter was set as a penalty currenlty set to 1.
2067 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2068 C second case jth atom is center
2069 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2070 C the third case kth atom is the center
2071 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2076 C write(iout,*)i,j,k,eij
2077 C The energy penalty calculated now time for the gradient part
2078 C derivative over rij
2079 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2080 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2085 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2086 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2089 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2090 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2092 C now derivative over rik
2093 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2094 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2099 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2100 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2103 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2104 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2106 C now derivative over rjk
2107 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2108 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2113 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2114 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2117 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2118 gvdwc(l,k)=gvdwc(l,k)+gg(l)