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'
103 include 'COMMON.SPLITELE'
111 double precision h_base
122 c integer itypi,itypj,k,l
123 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
124 double precision sig0ij,ljd,sig,fac,e1,e2
125 double precision dcosom1(3),dcosom2(3),ed
126 double precision pom1,pom2
127 double precision ljA,ljB,ljXs
128 double precision d_ljB(1:3)
129 double precision ssA,ssB,ssC,ssXs
130 double precision ssxm,ljxm,ssm,ljm
131 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
132 double precision f1,f2,h1,h2,hd1,hd2
133 double precision omega,delta_inv,deltasq_inv,fac1,fac2
135 double precision xm,d_xm(1:3)
136 double precision sslipi,sslipj,ssgradlipi,ssgradlipj
137 integer ici,icj,itypi,itypj
138 double precision boxshift,sscale,sscagrad
139 c-------END FIRST METHOD
140 c-------SECOND METHOD
141 c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
142 c-------END SECOND METHOD
145 logical checkstop,transgrad
146 common /sschecks/ checkstop,transgrad
148 integer icheck,nicheck,jcheck,njcheck
149 double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
150 c-------END TESTING CODE
157 c write (iout,*) "dyn_ssbond",resi,resj,ici,icj
159 if (ici.eq.0 .or. icj.eq.0) then
161 write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
162 & "Processor",me," attempt to create",
163 & " a disulfide link between non-cysteine residues ",restyp(i),i,
165 call MPI_Abort(MPI_COMM_WORLD,ierr)
167 write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
168 & "Processor",me," attempt to create",
169 & " a disulfide link between non-cysteine residues ",restyp(i),i,
175 dxi=dc_norm(1,nres+i)
176 dyi=dc_norm(2,nres+i)
177 dzi=dc_norm(3,nres+i)
178 dsci_inv=vbld_inv(i+nres)
182 call to_box(xi,yi,zi)
183 c write (iout,*) "After to_box i",xi,yi,zi
185 C define scaling factor for lipids
187 C if (positi.le.0) positi=positi+boxzsize
189 C first for peptide groups
190 c for each residue check if it is in lipid or lipid water border area
191 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
192 c write (iout,*) "After lipid_layer"
198 call to_box(xj,yj,zj)
199 c write (iout,*) "After to_box j",xj,yj,zj
201 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
202 c write (iout,*) "After lipid_layer"
204 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
205 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
206 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
207 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
208 xj=boxshift(xj-xi,boxxsize)
209 yj=boxshift(yj-yi,boxysize)
210 zj=boxshift(zj-zi,boxzsize)
211 c write (iout,*) "After boxshift"
213 dxj=dc_norm(1,nres+j)
214 dyj=dc_norm(2,nres+j)
215 dzj=dc_norm(3,nres+j)
216 dscj_inv=vbld_inv(j+nres)
218 chi1=chi(itypi,itypj)
219 chi2=chi(itypj,itypi)
226 alf12=0.5D0*(alf1+alf2)
228 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
229 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
230 sss=sscale((1.0d0/rij)/sigma(itypi,itypj),r_cut_int)
231 sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj),r_cut_int)
232 c The following are set in sc_angular
236 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
237 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
238 c om12=dxi*dxj+dyi*dyj+dzi*dzj
239 c write (iout,*) "Calling sc_angular"
242 c write (iout,*) "After sc_angular"
244 rij=1.0D0/rij ! Reset this so it makes sense
246 sig0ij=sigma(itypi,itypj)
247 sig=sig0ij*dsqrt(1.0D0/sigsq)
250 ljA=eps1*eps2rt**2*eps3rt**2
253 ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
258 deltat12=om2-om1+2.0d0
263 & +akth*(deltat1*deltat1+deltat2*deltat2)
264 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
265 ssxm=ssXs-0.5D0*ssB/ssA
268 c$$$c Some extra output
269 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
270 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
271 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
272 c$$$ if (ssx0.gt.0.0d0) then
273 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
277 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
278 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
279 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
281 c-------END TESTING CODE
284 c Stop and plot energy and derivative as a function of distance
285 c write (iout,*) "checkstop",checkstop
288 ssm=ssC-0.25D0*ssB*ssB/ssA
289 ljm=-0.25D0*ljB*bb/aa
291 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
299 if (.not.checkstop) then
306 if (checkstop) rij=(ssxm-1.0d0)+
307 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
308 c-------END TESTING CODE
310 if (rij.gt.ljxm) then
313 fac=(1.0D0/ljd)**expon
316 eij=eps1*eps2rt*eps3rt*(e1+e2)
319 eij=eij*eps2rt*eps3rt*sss
322 e1=e1*eps1*eps2rt**2*eps3rt**2
323 ed=-expon*(e1+eij)/ljd
325 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
326 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
327 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
328 eom12=eij*eps1_om12+eps2der*eps2rt_om12
329 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
330 else if (rij.lt.ssxm) then
333 eij=ssA*ssd*ssd+ssB*ssd+ssC
335 ed=2*akcm*ssd+akct*deltat12
336 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
338 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
339 eom1=-2*akth*deltat1-pom1-om2*pom2
340 eom2= 2*akth*deltat2+pom1-om1*pom2
343 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
345 d_ssxm(1)=0.5D0*akct/ssA
349 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
350 d_ljxm(2)=d_ljxm(1)*sigsq_om2
351 d_ljxm(3)=d_ljxm(1)*sigsq_om12
352 d_ljxm(1)=d_ljxm(1)*sigsq_om1
354 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
357 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
361 ssm=ssC-0.25D0*ssB*ssB/ssA
362 d_ssm(1)=0.5D0*akct*ssB/ssA
363 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
364 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
366 f1=(rij-xm)/(ssxm-xm)
367 f2=(rij-ssxm)/(xm-ssxm)
371 delta_inv=1.0d0/(xm-ssxm)
372 deltasq_inv=delta_inv*delta_inv
374 fac1=deltasq_inv*fac*(xm-rij)
375 fac2=deltasq_inv*fac*(rij-ssxm)
376 ed=delta_inv*(Ht*hd2-ssm*hd1)
378 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
379 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
380 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
381 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
384 ljm=-0.25D0*ljB*bb/aa
385 d_ljm(1)=-0.5D0*bb/aa*ljB
386 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
387 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
389 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
390 f1=(rij-ljxm)/(xm-ljxm)
391 f2=(rij-xm)/(ljxm-xm)
395 delta_inv=1.0d0/(ljxm-xm)
396 deltasq_inv=delta_inv*delta_inv
398 fac1=deltasq_inv*fac*(ljxm-rij)
399 fac2=deltasq_inv*fac*(rij-xm)
400 ed=delta_inv*(ljm*hd2-Ht*hd1)
402 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
403 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
404 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
405 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
407 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
409 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
415 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
416 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
417 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
419 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
420 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
421 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
422 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
425 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
427 c$$$ d_ljm(k)=ljm*d_ljB(k)
431 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
432 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
433 c$$$ d_ss(2)=akct*ssd
434 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
435 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
438 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
439 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
440 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
442 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
443 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
445 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
447 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
448 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
449 c$$$ h1=h_base(f1,hd1)
450 c$$$ h2=h_base(f2,hd2)
451 c$$$ eij=ss*h1+ljf*h2
452 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
453 c$$$ deltasq_inv=delta_inv*delta_inv
454 c$$$ fac=ljf*hd2-ss*hd1
455 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
456 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
457 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
458 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
459 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
460 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
461 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
463 c$$$ havebond=.false.
464 c$$$ if (ed.gt.0.0d0) havebond=.true.
465 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
469 c write (iout,*) "havebond",havebond
474 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
475 c write(iout,'(a15,f12.2,f8.1,2i5)')
476 c & "SSBOND_E_FORM",totT,t_bath,i,j
480 dyn_ssbond_ij(ici,icj)=eij
481 else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300)
483 dyn_ssbond_ij(ici,icj)=1.0d300
486 c write(iout,'(a15,f12.2,f8.1,2i5)')
487 c & "SSBOND_E_BREAK",totT,t_bath,i,j
494 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
495 & "CHECKSTOP",rij,eij,ed
500 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
507 c-------END TESTING CODE
508 gg_lipi(3)=ssgradlipi*eij
509 gg_lipj(3)=ssgradlipj*eij
512 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
513 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
516 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
519 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
520 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
521 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
522 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
523 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
524 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
528 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
533 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
534 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
540 C-----------------------------------------------------------------------------
542 double precision function h_base(x,deriv)
543 c A smooth function going 0->1 in range [0,1]
544 c It should NOT be called outside range [0,1], it will not work there.
551 double precision deriv
557 c Two parabolas put together. First derivative zero at extrema
558 c$$$ if (x.lt.0.5D0) then
559 c$$$ h_base=2.0D0*x*x
563 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
564 c$$$ deriv=4.0D0*deriv
567 c Third degree polynomial. First derivative zero at extrema
568 h_base=x*x*(3.0d0-2.0d0*x)
569 deriv=6.0d0*x*(1.0d0-x)
571 c Fifth degree polynomial. First and second derivatives zero at extrema
573 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
575 c$$$ deriv=deriv*deriv
576 c$$$ deriv=30.0d0*xsq*deriv
581 c----------------------------------------------------------------------------
583 subroutine dyn_set_nss
584 c Adjust nss and other relevant variables based on dyn_ssbond_ij
592 include 'COMMON.SBRIDGE'
593 include 'COMMON.CHAIN'
594 include 'COMMON.IOUNITS'
595 include 'COMMON.SETUP'
603 double precision emin
605 integer diff,allflag(maxss),allnss,
606 & allihpb(maxss),alljhpb(maxss),
607 & newnss,newihpb(maxss),newjhpb(maxss)
609 integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
610 integer g_newihpb(maxss),g_newjhpb(maxss),g_newnss
615 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
624 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
628 if (allflag(i).eq.0 .and.
629 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
630 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
634 if (emin.lt.1.0d300) then
637 if (allflag(i).eq.0 .and.
638 & (allihpb(i).eq.allihpb(imin) .or.
639 & alljhpb(i).eq.allihpb(imin) .or.
640 & allihpb(i).eq.alljhpb(imin) .or.
641 & alljhpb(i).eq.alljhpb(imin))) then
648 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
652 if (allflag(i).eq.1) then
654 newihpb(newnss)=allihpb(i)
655 newjhpb(newnss)=alljhpb(i)
660 if (nfgtasks.gt.1)then
662 call MPI_Reduce(newnss,g_newnss,1,
663 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
664 call MPI_Gather(newnss,1,MPI_INTEGER,
665 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
668 displ(i)=i_newnss(i-1)+displ(i-1)
670 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
671 & g_newihpb,i_newnss,displ,MPI_INTEGER,
673 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
674 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
676 if(fg_rank.eq.0) then
677 c print *,'g_newnss',g_newnss
678 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
679 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
682 newihpb(i)=g_newihpb(i)
683 newjhpb(i)=g_newjhpb(i)
691 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
696 if (idssb(i).eq.newihpb(j) .and.
697 & jdssb(i).eq.newjhpb(j)) found=.true.
701 c if (.not.found.and.fg_rank.eq.0)
702 c & write(iout,'(a15,f12.2,f8.1,2i5)')
703 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
711 if (newihpb(i).eq.idssb(j) .and.
712 & newjhpb(i).eq.jdssb(j)) found=.true.
716 c if (.not.found.and.fg_rank.eq.0)
717 c & write(iout,'(a15,f12.2,f8.1,2i5)')
718 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
733 C-----------------------------------------------------------------------------
734 C-----------------------------------------------------------------------------
735 C-----------------------------------------------------------------------------
737 c$$$c-----------------------------------------------------------------------------
739 c$$$ subroutine ss_relax(i_in,j_in)
743 c$$$ include 'DIMENSIONS'
744 c$$$ include 'COMMON.VAR'
745 c$$$ include 'COMMON.CHAIN'
746 c$$$ include 'COMMON.IOUNITS'
747 c$$$ include 'COMMON.INTERACT'
749 c$$$c Input arguments
750 c$$$ integer i_in,j_in
752 c$$$c Local variables
753 c$$$ integer i,iretcode,nfun_sc
755 c$$$ double precision var(maxvar),e_sc,etot
762 c$$$ mask_side(i_in)=1
763 c$$$ mask_side(j_in)=1
765 c$$$c Minimize the two selected side-chains
766 c$$$ call overlap_sc(scfail) ! Better not fail!
767 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
774 c$$$c-------------------------------------------------------------
776 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
777 c$$$c Minimize side-chains only, starting from geom but without modifying
779 c$$$c If mask_r is already set, only the selected side-chains are minimized,
780 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
784 c$$$ include 'DIMENSIONS'
785 c$$$ include 'COMMON.IOUNITS'
786 c$$$ include 'COMMON.VAR'
787 c$$$ include 'COMMON.CHAIN'
788 c$$$ include 'COMMON.GEO'
789 c$$$ include 'COMMON.MINIM'
791 c$$$ common /srutu/ icall
793 c$$$c Output arguments
794 c$$$ double precision etot_sc
795 c$$$ integer iretcode,nfun
797 c$$$c External functions/subroutines
798 c$$$ external func_sc,grad_sc,fdum
800 c$$$c Local variables
802 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
804 c$$$ double precision rdum(1)
805 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
807 c$$$ integer i,nvar_restr
810 c$$$cmc start_minim=.true.
811 c$$$ call deflt(2,iv,liv,lv,v)
812 c$$$* 12 means fresh start, dont call deflt
814 c$$$* max num of fun calls
815 c$$$ if (maxfun.eq.0) maxfun=500
817 c$$$* max num of iterations
818 c$$$ if (maxmin.eq.0) maxmin=1000
820 c$$$* controls output
822 c$$$* selects output unit
824 c$$$c iv(21)=iout ! DEBUG
825 c$$$c iv(21)=8 ! DEBUG
826 c$$$* 1 means to print out result
828 c$$$c iv(22)=1 ! DEBUG
829 c$$$* 1 means to print out summary stats
831 c$$$c iv(23)=1 ! DEBUG
832 c$$$* 1 means to print initial x and d
834 c$$$c iv(24)=1 ! DEBUG
835 c$$$* min val for v(radfac) default is 0.1
837 c$$$* max val for v(radfac) default is 4.0
840 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
841 c$$$* the sumsl default is 0.1
843 c$$$* false conv if (act fnctn decrease) .lt. v(34)
844 c$$$* the sumsl default is 100*machep
845 c$$$ v(34)=v(34)/100.0D0
846 c$$$* absolute convergence
847 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
849 c$$$* relative convergence
850 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
852 c$$$* controls initial step size
854 c$$$* large vals of d correspond to small components of step
858 c$$$ do i=nphi+1,nvar
862 c$$$ call geom_to_var(nvar,x)
863 c$$$ IF (mask_r) THEN
864 c$$$ do i=1,nres ! Just in case...
868 c$$$ call x2xx(x,xx,nvar_restr)
869 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
870 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
873 c$$$c When minimizing ALL side-chains, etotal_sc is a little
874 c$$$c faster if we don't set mask_r
880 c$$$ call x2xx(x,xx,nvar_restr)
881 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
882 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
885 c$$$ call var_to_geom(nvar,x)
886 c$$$ call chainbuild_sc
893 c$$$C--------------------------------------------------------------------------
895 c$$$ subroutine chainbuild_sc
897 c$$$ include 'DIMENSIONS'
898 c$$$ include 'COMMON.VAR'
899 c$$$ include 'COMMON.INTERACT'
901 c$$$c Local variables
906 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
907 c$$$ call locate_side_chain(i)
914 c$$$C--------------------------------------------------------------------------
916 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
920 c$$$ include 'DIMENSIONS'
921 c$$$ include 'COMMON.DERIV'
922 c$$$ include 'COMMON.VAR'
923 c$$$ include 'COMMON.MINIM'
924 c$$$ include 'COMMON.IOUNITS'
926 c$$$c Input arguments
928 c$$$ double precision x(maxvar)
929 c$$$ double precision ufparm
932 c$$$c Input/Output arguments
934 c$$$ integer uiparm(1)
935 c$$$ double precision urparm(1)
937 c$$$c Output arguments
938 c$$$ double precision f
940 c$$$c Local variables
941 c$$$ double precision energia(0:n_ene)
943 c$$$c Variables used to intercept NaNs
944 c$$$ double precision x_sum
953 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
956 c$$$ x_sum=x_sum+x(i_NAN)
958 c$$$c Calculate the energy only if the coordinates are ok
959 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
960 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
966 c$$$ call var_to_geom_restr(n,x)
968 c$$$ call chainbuild_sc
969 c$$$ call etotal_sc(energia(0))
971 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
980 c$$$c-------------------------------------------------------
982 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
986 c$$$ include 'DIMENSIONS'
987 c$$$ include 'COMMON.CHAIN'
988 c$$$ include 'COMMON.DERIV'
989 c$$$ include 'COMMON.VAR'
990 c$$$ include 'COMMON.INTERACT'
991 c$$$ include 'COMMON.MINIM'
993 c$$$c Input arguments
995 c$$$ double precision x(maxvar)
996 c$$$ double precision ufparm
999 c$$$c Input/Output arguments
1001 c$$$ integer uiparm(1)
1002 c$$$ double precision urparm(1)
1004 c$$$c Output arguments
1005 c$$$ double precision g(maxvar)
1007 c$$$c Local variables
1008 c$$$ double precision f,gphii,gthetai,galphai,gomegai
1009 c$$$ integer ig,ind,i,j,k,igall,ij
1012 c$$$ icg=mod(nf,2)+1
1013 c$$$ if (nf-nfl+1) 20,30,40
1014 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1015 c$$$c write (iout,*) 'grad 20'
1016 c$$$ if (nf.eq.0) return
1018 c$$$ 30 call var_to_geom_restr(n,x)
1019 c$$$ call chainbuild_sc
1021 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1023 c$$$ 40 call cartder
1025 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1031 c$$$ IF (mask_phi(i+2).eq.1) THEN
1033 c$$$ do j=i+1,nres-1
1036 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1037 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1043 c$$$ ind=ind+nres-1-i
1050 c$$$ IF (mask_theta(i+2).eq.1) THEN
1053 c$$$ do j=i+1,nres-1
1056 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1057 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1062 c$$$ ind=ind+nres-1-i
1067 c$$$ if (itype(i).ne.10) then
1068 c$$$ IF (mask_side(i).eq.1) THEN
1072 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1081 c$$$ if (itype(i).ne.10) then
1082 c$$$ IF (mask_side(i).eq.1) THEN
1086 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1094 c$$$C Add the components corresponding to local energy terms.
1101 c$$$ if (mask_phi(i).eq.1) then
1103 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1109 c$$$ if (mask_theta(i).eq.1) then
1111 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1117 c$$$ if (itype(i).ne.10) then
1119 c$$$ if (mask_side(i).eq.1) then
1121 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1128 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1134 c$$$C-----------------------------------------------------------------------------
1136 c$$$ subroutine etotal_sc(energy_sc)
1140 c$$$ include 'DIMENSIONS'
1141 c$$$ include 'COMMON.VAR'
1142 c$$$ include 'COMMON.INTERACT'
1143 c$$$ include 'COMMON.DERIV'
1144 c$$$ include 'COMMON.FFIELD'
1146 c$$$c Output arguments
1147 c$$$ double precision energy_sc(0:n_ene)
1149 c$$$c Local variables
1150 c$$$ double precision evdw,escloc
1155 c$$$ energy_sc(i)=0.0D0
1158 c$$$ if (mask_r) then
1159 c$$$ call egb_sc(evdw)
1160 c$$$ call esc_sc(escloc)
1163 c$$$ call esc(escloc)
1166 c$$$ if (evdw.eq.1.0D20) then
1167 c$$$ energy_sc(0)=evdw
1169 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1171 c$$$ energy_sc(1)=evdw
1172 c$$$ energy_sc(12)=escloc
1175 c$$$C Sum up the components of the Cartesian gradient.
1179 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1186 c$$$C-----------------------------------------------------------------------------
1188 c$$$ subroutine egb_sc(evdw)
1190 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1191 c$$$C assuming the Gay-Berne potential of interaction.
1193 c$$$ implicit real*8 (a-h,o-z)
1194 c$$$ include 'DIMENSIONS'
1195 c$$$ include 'COMMON.GEO'
1196 c$$$ include 'COMMON.VAR'
1197 c$$$ include 'COMMON.LOCAL'
1198 c$$$ include 'COMMON.CHAIN'
1199 c$$$ include 'COMMON.DERIV'
1200 c$$$ include 'COMMON.NAMES'
1201 c$$$ include 'COMMON.INTERACT'
1202 c$$$ include 'COMMON.IOUNITS'
1203 c$$$ include 'COMMON.CALC'
1204 c$$$ include 'COMMON.CONTROL'
1207 c$$$ energy_dec=.false.
1208 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1211 c$$$c if (icall.eq.0) lprn=.false.
1213 c$$$ do i=iatsc_s,iatsc_e
1215 c$$$ itypi1=itype(i+1)
1219 c$$$ dxi=dc_norm(1,nres+i)
1220 c$$$ dyi=dc_norm(2,nres+i)
1221 c$$$ dzi=dc_norm(3,nres+i)
1222 c$$$c dsci_inv=dsc_inv(itypi)
1223 c$$$ dsci_inv=vbld_inv(i+nres)
1224 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1225 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1227 c$$$C Calculate SC interaction energy.
1229 c$$$ do iint=1,nint_gr(i)
1230 c$$$ do j=istart(i,iint),iend(i,iint)
1231 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1234 c$$$c dscj_inv=dsc_inv(itypj)
1235 c$$$ dscj_inv=vbld_inv(j+nres)
1236 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1237 c$$$c & 1.0d0/vbld(j+nres)
1238 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1239 c$$$ sig0ij=sigma(itypi,itypj)
1240 c$$$ chi1=chi(itypi,itypj)
1241 c$$$ chi2=chi(itypj,itypi)
1242 c$$$ chi12=chi1*chi2
1243 c$$$ chip1=chip(itypi)
1244 c$$$ chip2=chip(itypj)
1245 c$$$ chip12=chip1*chip2
1246 c$$$ alf1=alp(itypi)
1247 c$$$ alf2=alp(itypj)
1248 c$$$ alf12=0.5D0*(alf1+alf2)
1249 c$$$C For diagnostics only!!!
1259 c$$$ xj=c(1,nres+j)-xi
1260 c$$$ yj=c(2,nres+j)-yi
1261 c$$$ zj=c(3,nres+j)-zi
1262 c$$$ dxj=dc_norm(1,nres+j)
1263 c$$$ dyj=dc_norm(2,nres+j)
1264 c$$$ dzj=dc_norm(3,nres+j)
1265 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1266 c$$$c write (iout,*) "j",j," dc_norm",
1267 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1268 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1269 c$$$ rij=dsqrt(rrij)
1270 c$$$C Calculate angle-dependent terms of energy and contributions to their
1272 c$$$ call sc_angular
1273 c$$$ sigsq=1.0D0/sigsq
1274 c$$$ sig=sig0ij*dsqrt(sigsq)
1275 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1276 c$$$c for diagnostics; uncomment
1277 c$$$c rij_shift=1.2*sig0ij
1278 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1279 c$$$ if (rij_shift.le.0.0D0) then
1281 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1282 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1283 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1286 c$$$ sigder=-sig*sigsq
1287 c$$$c---------------------------------------------------------------
1288 c$$$ rij_shift=1.0D0/rij_shift
1289 c$$$ fac=rij_shift**expon
1290 c$$$ e1=fac*fac*aa(itypi,itypj)
1291 c$$$ e2=fac*bb(itypi,itypj)
1292 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1293 c$$$ eps2der=evdwij*eps3rt
1294 c$$$ eps3der=evdwij*eps2rt
1295 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1296 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1297 c$$$ evdwij=evdwij*eps2rt*eps3rt
1298 c$$$ evdw=evdw+evdwij
1300 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1301 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1302 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1303 c$$$ & restyp(itypi),i,restyp(itypj),j,
1304 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1305 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1306 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1310 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1311 c$$$ & 'evdw',i,j,evdwij
1313 c$$$C Calculate gradient components.
1314 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1315 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1316 c$$$ sigder=fac*sigder
1319 c$$$C Calculate the radial part of the gradient
1323 c$$$C Calculate angular part of the gradient.
1329 c$$$ energy_dec=.false.
1333 c$$$c-----------------------------------------------------------------------------
1335 c$$$ subroutine esc_sc(escloc)
1336 c$$$C Calculate the local energy of a side chain and its derivatives in the
1337 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1338 c$$$C ALPHA and OMEGA.
1339 c$$$ implicit real*8 (a-h,o-z)
1340 c$$$ include 'DIMENSIONS'
1341 c$$$ include 'COMMON.GEO'
1342 c$$$ include 'COMMON.LOCAL'
1343 c$$$ include 'COMMON.VAR'
1344 c$$$ include 'COMMON.INTERACT'
1345 c$$$ include 'COMMON.DERIV'
1346 c$$$ include 'COMMON.CHAIN'
1347 c$$$ include 'COMMON.IOUNITS'
1348 c$$$ include 'COMMON.NAMES'
1349 c$$$ include 'COMMON.FFIELD'
1350 c$$$ include 'COMMON.CONTROL'
1351 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1352 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1353 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1354 c$$$ delta=0.02d0*pi
1356 c$$$c write (iout,'(a)') 'ESC'
1357 c$$$ do i=loc_start,loc_end
1358 c$$$ IF (mask_side(i).eq.1) THEN
1360 c$$$ if (it.eq.10) goto 1
1361 c$$$ nlobit=nlob(it)
1362 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1363 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1364 c$$$ theti=theta(i+1)-pipol
1365 c$$$ x(1)=dtan(theti)
1369 c$$$ if (x(2).gt.pi-delta) then
1371 c$$$ xtemp(2)=pi-delta
1373 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1375 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1376 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1377 c$$$ & escloci,dersc(2))
1378 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1379 c$$$ & ddersc0(1),dersc(1))
1380 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1381 c$$$ & ddersc0(3),dersc(3))
1382 c$$$ xtemp(2)=pi-delta
1383 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1385 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1386 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1387 c$$$ & dersc0(2),esclocbi,dersc02)
1388 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1389 c$$$ & dersc12,dersc01)
1390 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1391 c$$$ dersc0(1)=dersc01
1392 c$$$ dersc0(2)=dersc02
1393 c$$$ dersc0(3)=0.0d0
1395 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1397 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1398 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1399 c$$$c & esclocbi,ss,ssd
1400 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1401 c$$$c escloci=esclocbi
1402 c$$$c write (iout,*) escloci
1403 c$$$ else if (x(2).lt.delta) then
1407 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1409 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1410 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1411 c$$$ & escloci,dersc(2))
1412 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1413 c$$$ & ddersc0(1),dersc(1))
1414 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1415 c$$$ & ddersc0(3),dersc(3))
1417 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1419 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1420 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1421 c$$$ & dersc0(2),esclocbi,dersc02)
1422 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1423 c$$$ & dersc12,dersc01)
1424 c$$$ dersc0(1)=dersc01
1425 c$$$ dersc0(2)=dersc02
1426 c$$$ dersc0(3)=0.0d0
1427 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1429 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1431 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1432 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1433 c$$$c & esclocbi,ss,ssd
1434 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1435 c$$$c write (iout,*) escloci
1437 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1440 c$$$ escloc=escloc+escloci
1441 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1442 c$$$ & 'escloc',i,escloci
1443 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1445 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1446 c$$$ & wscloc*dersc(1)
1447 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1448 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1455 c$$$C-----------------------------------------------------------------------------
1457 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1459 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1460 c$$$C assuming the Gay-Berne potential of interaction.
1462 c$$$ implicit real*8 (a-h,o-z)
1463 c$$$ include 'DIMENSIONS'
1464 c$$$ include 'COMMON.GEO'
1465 c$$$ include 'COMMON.VAR'
1466 c$$$ include 'COMMON.LOCAL'
1467 c$$$ include 'COMMON.CHAIN'
1468 c$$$ include 'COMMON.DERIV'
1469 c$$$ include 'COMMON.NAMES'
1470 c$$$ include 'COMMON.INTERACT'
1471 c$$$ include 'COMMON.IOUNITS'
1472 c$$$ include 'COMMON.CALC'
1473 c$$$ include 'COMMON.CONTROL'
1476 c$$$ energy_dec=.false.
1477 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1481 c$$$c$$$ do i=iatsc_s,iatsc_e
1484 c$$$ itypi1=itype(i+1)
1488 c$$$ dxi=dc_norm(1,nres+i)
1489 c$$$ dyi=dc_norm(2,nres+i)
1490 c$$$ dzi=dc_norm(3,nres+i)
1491 c$$$c dsci_inv=dsc_inv(itypi)
1492 c$$$ dsci_inv=vbld_inv(i+nres)
1493 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1494 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1496 c$$$C Calculate SC interaction energy.
1498 c$$$c$$$ do iint=1,nint_gr(i)
1499 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1503 c$$$c dscj_inv=dsc_inv(itypj)
1504 c$$$ dscj_inv=vbld_inv(j+nres)
1505 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1506 c$$$c & 1.0d0/vbld(j+nres)
1507 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1508 c$$$ sig0ij=sigma(itypi,itypj)
1509 c$$$ chi1=chi(itypi,itypj)
1510 c$$$ chi2=chi(itypj,itypi)
1511 c$$$ chi12=chi1*chi2
1512 c$$$ chip1=chip(itypi)
1513 c$$$ chip2=chip(itypj)
1514 c$$$ chip12=chip1*chip2
1515 c$$$ alf1=alp(itypi)
1516 c$$$ alf2=alp(itypj)
1517 c$$$ alf12=0.5D0*(alf1+alf2)
1518 c$$$C For diagnostics only!!!
1528 c$$$ xj=c(1,nres+j)-xi
1529 c$$$ yj=c(2,nres+j)-yi
1530 c$$$ zj=c(3,nres+j)-zi
1531 c$$$ dxj=dc_norm(1,nres+j)
1532 c$$$ dyj=dc_norm(2,nres+j)
1533 c$$$ dzj=dc_norm(3,nres+j)
1534 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1535 c$$$c write (iout,*) "j",j," dc_norm",
1536 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1537 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1538 c$$$ rij=dsqrt(rrij)
1539 c$$$C Calculate angle-dependent terms of energy and contributions to their
1541 c$$$ call sc_angular
1542 c$$$ sigsq=1.0D0/sigsq
1543 c$$$ sig=sig0ij*dsqrt(sigsq)
1544 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1545 c$$$c for diagnostics; uncomment
1546 c$$$c rij_shift=1.2*sig0ij
1547 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1548 c$$$ if (rij_shift.le.0.0D0) then
1550 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1551 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1552 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1555 c$$$ sigder=-sig*sigsq
1556 c$$$c---------------------------------------------------------------
1557 c$$$ rij_shift=1.0D0/rij_shift
1558 c$$$ fac=rij_shift**expon
1559 c$$$ e1=fac*fac*aa(itypi,itypj)
1560 c$$$ e2=fac*bb(itypi,itypj)
1561 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1562 c$$$ eps2der=evdwij*eps3rt
1563 c$$$ eps3der=evdwij*eps2rt
1564 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1565 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1566 c$$$ evdwij=evdwij*eps2rt*eps3rt
1567 c$$$ evdw=evdw+evdwij
1569 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1570 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1571 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1572 c$$$ & restyp(itypi),i,restyp(itypj),j,
1573 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1574 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1575 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1579 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1580 c$$$ & 'evdw',i,j,evdwij
1582 c$$$C Calculate gradient components.
1583 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1584 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1585 c$$$ sigder=fac*sigder
1588 c$$$C Calculate the radial part of the gradient
1592 c$$$C Calculate angular part of the gradient.
1595 c$$$c$$$ enddo ! iint
1597 c$$$ energy_dec=.false.
1601 c$$$C-----------------------------------------------------------------------------
1603 c$$$ subroutine perturb_side_chain(i,angle)
1607 c$$$ include 'DIMENSIONS'
1608 c$$$ include 'COMMON.CHAIN'
1609 c$$$ include 'COMMON.GEO'
1610 c$$$ include 'COMMON.VAR'
1611 c$$$ include 'COMMON.LOCAL'
1612 c$$$ include 'COMMON.IOUNITS'
1614 c$$$c External functions
1615 c$$$ external ran_number
1616 c$$$ double precision ran_number
1618 c$$$c Input arguments
1620 c$$$ double precision angle ! In degrees
1622 c$$$c Local variables
1624 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1628 c$$$ rad_ang=angle*deg2rad
1631 c$$$ do while (length.lt.0.01)
1632 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1633 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1634 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1635 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1636 c$$$ + rand_v(3)*rand_v(3)
1637 c$$$ length=sqrt(length)
1638 c$$$ rand_v(1)=rand_v(1)/length
1639 c$$$ rand_v(2)=rand_v(2)/length
1640 c$$$ rand_v(3)=rand_v(3)/length
1641 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1642 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1643 c$$$ length=1.0D0-cost*cost
1644 c$$$ if (length.lt.0.0D0) length=0.0D0
1645 c$$$ length=sqrt(length)
1646 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1647 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1648 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1650 c$$$ rand_v(1)=rand_v(1)/length
1651 c$$$ rand_v(2)=rand_v(2)/length
1652 c$$$ rand_v(3)=rand_v(3)/length
1654 c$$$ cost=dcos(rad_ang)
1655 c$$$ sint=dsin(rad_ang)
1656 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1657 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1658 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1659 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1660 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1661 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1662 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1663 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1664 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1666 c$$$ call chainbuild_cart
1671 c$$$c----------------------------------------------------------------------------
1673 c$$$ subroutine ss_relax3(i_in,j_in)
1677 c$$$ include 'DIMENSIONS'
1678 c$$$ include 'COMMON.VAR'
1679 c$$$ include 'COMMON.CHAIN'
1680 c$$$ include 'COMMON.IOUNITS'
1681 c$$$ include 'COMMON.INTERACT'
1683 c$$$c External functions
1684 c$$$ external ran_number
1685 c$$$ double precision ran_number
1687 c$$$c Input arguments
1688 c$$$ integer i_in,j_in
1690 c$$$c Local variables
1691 c$$$ double precision energy_sc(0:n_ene),etot
1692 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1693 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1694 c$$$ integer n,i_pert,i
1695 c$$$ logical notdone
1704 c$$$ mask_side(i_in)=1
1705 c$$$ mask_side(j_in)=1
1707 c$$$ call etotal_sc(energy_sc)
1708 c$$$ etot=energy_sc(0)
1709 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1710 c$$$c + energy_sc(1),energy_sc(12)
1714 c$$$ do while (notdone)
1715 c$$$ if (mod(n,2).eq.0) then
1723 c$$$ org_dc(i)=dc(i,i_pert+nres)
1724 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1725 c$$$ org_c(i)=c(i,i_pert+nres)
1727 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1728 c$$$ call perturb_side_chain(i_pert,ang_pert)
1729 c$$$ call etotal_sc(energy_sc)
1730 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1731 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1732 c$$$ if (rand_fact.lt.exp_fact) then
1733 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1734 c$$$c + energy_sc(1),energy_sc(12)
1735 c$$$ etot=energy_sc(0)
1737 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1738 c$$$c + energy_sc(1),energy_sc(12)
1740 c$$$ dc(i,i_pert+nres)=org_dc(i)
1741 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1742 c$$$ c(i,i_pert+nres)=org_c(i)
1746 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1754 c$$$c----------------------------------------------------------------------------
1756 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1758 c$$$ include 'DIMENSIONS'
1760 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1761 c$$$*********************************************************************
1762 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1763 c$$$* the calling subprogram. *
1764 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1765 c$$$* calculated in the usual pythagorean way. *
1766 c$$$* absolute convergence occurs when the function is within v(31) of *
1767 c$$$* zero. unless you know the minimum value in advance, abs convg *
1768 c$$$* is probably not useful. *
1769 c$$$* relative convergence is when the model predicts that the function *
1770 c$$$* will decrease by less than v(32)*abs(fun). *
1771 c$$$*********************************************************************
1772 c$$$ include 'COMMON.IOUNITS'
1773 c$$$ include 'COMMON.VAR'
1774 c$$$ include 'COMMON.GEO'
1775 c$$$ include 'COMMON.MINIM'
1776 c$$$ include 'COMMON.CHAIN'
1778 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1779 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1780 c$$$ + orig_ss_dist(maxres2,maxres2)
1782 c$$$ double precision etot
1783 c$$$ integer iretcode,nfun,i_in,j_in
1786 c$$$ double precision dist
1787 c$$$ external ss_func,fdum
1788 c$$$ double precision ss_func,fdum
1790 c$$$ integer iv(liv),uiparm(2)
1791 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1795 c$$$ call deflt(2,iv,liv,lv,v)
1796 c$$$* 12 means fresh start, dont call deflt
1798 c$$$* max num of fun calls
1799 c$$$ if (maxfun.eq.0) maxfun=500
1801 c$$$* max num of iterations
1802 c$$$ if (maxmin.eq.0) maxmin=1000
1804 c$$$* controls output
1806 c$$$* selects output unit
1809 c$$$* 1 means to print out result
1811 c$$$* 1 means to print out summary stats
1813 c$$$* 1 means to print initial x and d
1815 c$$$* min val for v(radfac) default is 0.1
1817 c$$$* max val for v(radfac) default is 4.0
1820 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1821 c$$$* the sumsl default is 0.1
1823 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1824 c$$$* the sumsl default is 100*machep
1825 c$$$ v(34)=v(34)/100.0D0
1826 c$$$* absolute convergence
1827 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1830 c$$$* relative convergence
1831 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1834 c$$$* controls initial step size
1836 c$$$* large vals of d correspond to small components of step
1843 c$$$ orig_ss_dc(j,i)=dc(j,i)
1846 c$$$ call geom_to_var(nvar,orig_ss_var)
1850 c$$$ orig_ss_dist(j,i)=dist(j,i)
1851 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1852 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1853 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1865 c$$$ if (ialph(i,1).gt.0) then
1868 c$$$ x(k)=dc(j,i+nres)
1875 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1878 c$$$ nfun=iv(6)+iv(30)
1888 c$$$ if (ialph(i,1).gt.0) then
1891 c$$$ dc(j,i+nres)=x(k)
1895 c$$$ call chainbuild_cart
1900 c$$$C-----------------------------------------------------------------------------
1902 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1904 c$$$ include 'DIMENSIONS'
1905 c$$$ include 'COMMON.DERIV'
1906 c$$$ include 'COMMON.IOUNITS'
1907 c$$$ include 'COMMON.VAR'
1908 c$$$ include 'COMMON.CHAIN'
1909 c$$$ include 'COMMON.INTERACT'
1910 c$$$ include 'COMMON.SBRIDGE'
1912 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1913 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1914 c$$$ + orig_ss_dist(maxres2,maxres2)
1917 c$$$ double precision x(maxres6)
1919 c$$$ double precision f
1920 c$$$ integer uiparm(2)
1921 c$$$ real*8 urparm(1)
1922 c$$$ external ufparm
1923 c$$$ double precision ufparm
1926 c$$$ double precision dist
1928 c$$$ integer i,j,k,ss_i,ss_j
1929 c$$$ double precision tempf,var(maxvar)
1944 c$$$ if (ialph(i,1).gt.0) then
1947 c$$$ dc(j,i+nres)=x(k)
1951 c$$$ call chainbuild_cart
1953 c$$$ call geom_to_var(nvar,var)
1955 c$$$c Constraints on all angles
1957 c$$$ tempf=var(i)-orig_ss_var(i)
1958 c$$$ f=f+tempf*tempf
1961 c$$$c Constraints on all distances
1963 c$$$ if (i.gt.1) then
1964 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1965 c$$$ f=f+tempf*tempf
1968 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
1969 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
1970 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
1971 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1972 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
1973 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1974 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
1975 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
1979 c$$$c Constraints for the relevant CYS-CYS
1980 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
1981 c$$$ f=f+tempf*tempf
1982 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
1984 c$$$c$$$ if (nf.ne.nfl) then
1985 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
1986 c$$$c$$$ + f,dist(5+nres,14+nres)
1994 c$$$C-----------------------------------------------------------------------------
1995 c$$$C-----------------------------------------------------------------------------
1996 subroutine triple_ssbond_ene(resi,resj,resk,eij)
1997 include 'DIMENSIONS'
1998 include 'COMMON.SBRIDGE'
1999 include 'COMMON.CHAIN'
2000 include 'COMMON.DERIV'
2001 include 'COMMON.LOCAL'
2002 include 'COMMON.INTERACT'
2003 include 'COMMON.VAR'
2004 include 'COMMON.IOUNITS'
2005 include 'COMMON.CALC'
2012 c External functions
2013 double precision h_base
2017 integer resi,resj,resk
2020 double precision eij,eij1,eij2,eij3
2024 c integer itypi,itypj,k,l
2025 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2026 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2027 double precision xik,yik,zik,xjk,yjk,zjk
2028 double precision sig0ij,ljd,sig,fac,e1,e2
2029 double precision dcosom1(3),dcosom2(3),ed
2030 double precision pom1,pom2
2031 double precision ljA,ljB,ljXs
2032 double precision d_ljB(1:3)
2033 double precision ssA,ssB,ssC,ssXs
2034 double precision ssxm,ljxm,ssm,ljm
2035 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2036 if (dtriss.eq.0) return
2040 C write(iout,*) resi,resj,resk
2042 dxi=dc_norm(1,nres+i)
2043 dyi=dc_norm(2,nres+i)
2044 dzi=dc_norm(3,nres+i)
2045 dsci_inv=vbld_inv(i+nres)
2055 dxj=dc_norm(1,nres+j)
2056 dyj=dc_norm(2,nres+j)
2057 dzj=dc_norm(3,nres+j)
2058 dscj_inv=vbld_inv(j+nres)
2064 dxk=dc_norm(1,nres+k)
2065 dyk=dc_norm(2,nres+k)
2066 dzk=dc_norm(3,nres+k)
2067 dscj_inv=vbld_inv(k+nres)
2077 rrij=(xij*xij+yij*yij+zij*zij)
2078 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2079 rrik=(xik*xik+yik*yik+zik*zik)
2081 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2083 C there are three combination of distances for each trisulfide bonds
2084 C The first case the ith atom is the center
2085 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2086 C distance y is second distance the a,b,c,d are parameters derived for
2087 C this problem d parameter was set as a penalty currenlty set to 1.
2088 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2089 C second case jth atom is center
2090 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2091 C the third case kth atom is the center
2092 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2097 C write(iout,*)i,j,k,eij
2098 C The energy penalty calculated now time for the gradient part
2099 C derivative over rij
2100 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2101 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2106 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2107 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2110 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2111 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2113 C now derivative over rik
2114 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2115 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2120 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2121 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2124 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2125 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2127 C now derivative over rjk
2128 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2129 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2134 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2135 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2138 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2139 gvdwc(l,k)=gvdwc(l,k)+gg(l)