1 c----------------------------------------------------------------------------
2 subroutine check_energies
9 include 'COMMON.IOUNITS'
10 include 'COMMON.SBRIDGE'
11 include 'COMMON.LOCAL'
15 double precision ran_number
19 integer i,j,k,l,lmax,p,pmax
20 double precision rmin,rmax
24 double precision wi,rij,tj,pj
48 ct wi=ran_number(0.0D0,pi)
49 c wi=ran_number(0.0D0,pi/6.0D0)
51 ct tj=ran_number(0.0D0,pi)
52 ct pj=ran_number(0.0D0,pi)
53 c pj=ran_number(0.0D0,pi/6.0D0)
57 ct rij=ran_number(rmin,rmax)
59 c(1,j)=d*sin(pj)*cos(tj)
60 c(2,j)=d*sin(pj)*sin(tj)
69 dc(k,nres+i)=c(k,nres+i)-c(k,i)
70 dc_norm(k,nres+i)=dc(k,nres+i)/d
71 dc(k,nres+j)=c(k,nres+j)-c(k,j)
72 dc_norm(k,nres+j)=dc(k,nres+j)/d
75 call dyn_ssbond_ene(i,j,eij)
84 C-----------------------------------------------------------------------------
86 subroutine dyn_ssbond_ene(resi,resj,eij)
91 include 'COMMON.SBRIDGE'
92 include 'COMMON.CHAIN'
93 include 'COMMON.DERIV'
94 include 'COMMON.LOCAL'
95 include 'COMMON.INTERACT'
97 include 'COMMON.IOUNITS'
106 double precision h_base
117 c integer itypi,itypj,k,l
118 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
119 double precision sig0ij,ljd,sig,fac,e1,e2
120 double precision dcosom1(3),dcosom2(3),ed
121 double precision pom1,pom2
122 double precision ljA,ljB,ljXs
123 double precision d_ljB(1:3)
124 double precision ssA,ssB,ssC,ssXs
125 double precision ssxm,ljxm,ssm,ljm
126 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
127 double precision f1,f2,h1,h2,hd1,hd2
128 double precision omega,delta_inv,deltasq_inv,fac1,fac2
130 double precision xm,d_xm(1:3)
131 c-------END FIRST METHOD
132 c-------SECOND METHOD
133 c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
134 c-------END SECOND METHOD
137 logical checkstop,transgrad
138 common /sschecks/ checkstop,transgrad
140 integer icheck,nicheck,jcheck,njcheck
141 double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
142 c-------END TESTING CODE
149 dxi=dc_norm(1,nres+i)
150 dyi=dc_norm(2,nres+i)
151 dzi=dc_norm(3,nres+i)
152 dsci_inv=vbld_inv(i+nres)
157 if (xi.lt.0) xi=xi+boxxsize
159 if (yi.lt.0) yi=yi+boxysize
161 if (zi.lt.0) zi=zi+boxzsize
162 C define scaling factor for lipids
164 C if (positi.le.0) positi=positi+boxzsize
166 C first for peptide groups
167 c for each residue check if it is in lipid or lipid water border area
168 if ((zi.gt.bordlipbot)
169 &.and.(zi.lt.bordliptop)) then
170 C the energy transfer exist
171 if (zi.lt.buflipbot) then
172 C what fraction I am in
174 & ((positi-bordlipbot)/lipbufthick)
175 C lipbufthick is thickenes of lipid buffore
176 sslipi=sscalelip(fracinbuf)
177 ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
178 elseif (zi.gt.bufliptop) then
179 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
180 sslipi=sscalelip(fracinbuf)
181 ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
195 if (xj.lt.0) xj=xj+boxxsize
197 if (yj.lt.0) yj=yj+boxysize
199 if (zj.lt.0) zj=zj+boxzsize
200 if ((zj.gt.bordlipbot)
201 &.and.(zj.lt.bordliptop)) then
202 C the energy transfer exist
203 if (zj.lt.buflipbot) then
204 C what fraction I am in
206 & ((positi-bordlipbot)/lipbufthick)
207 C lipbufthick is thickenes of lipid buffore
208 sslipj=sscalelip(fracinbuf)
209 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
210 elseif (zi.gt.bufliptop) then
211 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
212 sslipj=sscalelip(fracinbuf)
213 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
222 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
223 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
224 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
225 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
227 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
228 c write (iout,*) "init",xi,yi,xi,xj,yj,zj,sqrt(dist_init)
240 xj=xj_safe+xshift*boxxsize
241 yj=yj_safe+yshift*boxysize
242 zj=zj_safe+zshift*boxzsize
243 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
244 c write (iout,*) xshift,yshift,zshift,
245 c & " dist_temp",sqrt(dist_temp)
246 if(dist_temp.lt.dist_init) then
256 c write (iout,*) "subchap",subchap,xj_temp,yj_temp,zj_temp,
257 c & xj_safe,yj_safe,zj_safe
258 if (subchap.eq.1) then
267 c write (iout,*) "diff",xj,yj,zj
269 C xj=c(1,nres+j)-c(1,nres+i)
270 C yj=c(2,nres+j)-c(2,nres+i)
271 C zj=c(3,nres+j)-c(3,nres+i)
272 dxj=dc_norm(1,nres+j)
273 dyj=dc_norm(2,nres+j)
274 dzj=dc_norm(3,nres+j)
275 dscj_inv=vbld_inv(j+nres)
277 chi1=chi(itypi,itypj)
278 chi2=chi(itypj,itypi)
285 alf12=0.5D0*(alf1+alf2)
287 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
288 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
289 c sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
290 c sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
293 c write (iout,*) xj,yj,zj,1.0d0/rij
294 c The following are set in sc_angular
298 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
299 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
300 c om12=dxi*dxj+dyi*dyj+dzi*dzj
302 rij=1.0D0/rij ! Reset this so it makes sense
304 sig0ij=sigma(itypi,itypj)
305 sig=sig0ij*dsqrt(1.0D0/sigsq)
308 ljA=eps1*eps2rt**2*eps3rt**2
311 ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
316 deltat12=om2-om1+2.0d0
321 & +akth*(deltat1*deltat1+deltat2*deltat2)
322 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
323 ssxm=ssXs-0.5D0*ssB/ssA
324 c write (iout,*)"akcm",akcm," akct",akct," delta12",delta12,
325 c & " ss_depth",ss_depth," v1ss",v1ss," v2ss",v2ss," v3ss",v3ss
328 c$$$c Some extra output
329 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
330 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
331 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
332 c$$$ if (ssx0.gt.0.0d0) then
333 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
337 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
338 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
339 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
341 c-------END TESTING CODE
344 c Stop and plot energy and derivative as a function of distance
346 ssm=ssC-0.25D0*ssB*ssB/ssA
347 ljm=-0.25D0*ljB*bb/aa
349 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
357 if (.not.checkstop) then
364 if (checkstop) rij=(ssxm-1.0d0)+
365 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
366 c-------END TESTING CODE
368 if (rij.gt.ljxm) then
371 fac=(1.0D0/ljd)**expon
374 eij=eps1*eps2rt*eps3rt*(e1+e2)
377 eij=eij*eps2rt*eps3rt*sss
380 e1=e1*eps1*eps2rt**2*eps3rt**2
381 ed=-expon*(e1+eij)/ljd
383 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
384 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
385 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
386 eom12=eij*eps1_om12+eps2der*eps2rt_om12
387 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
388 else if (rij.lt.ssxm) then
391 eij=ssA*ssd*ssd+ssB*ssd+ssC
393 ed=2*akcm*ssd+akct*deltat12
394 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
396 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
397 eom1=-2*akth*deltat1-pom1-om2*pom2
398 eom2= 2*akth*deltat2+pom1-om1*pom2
401 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
403 d_ssxm(1)=0.5D0*akct/ssA
407 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
408 d_ljxm(2)=d_ljxm(1)*sigsq_om2
409 d_ljxm(3)=d_ljxm(1)*sigsq_om12
410 d_ljxm(1)=d_ljxm(1)*sigsq_om1
412 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
415 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
419 ssm=ssC-0.25D0*ssB*ssB/ssA
420 c write(iout,*) "ssC",ssC," ssB",ssB," ssA",ssA
421 d_ssm(1)=0.5D0*akct*ssB/ssA
422 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
423 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
425 f1=(rij-xm)/(ssxm-xm)
426 f2=(rij-ssxm)/(xm-ssxm)
429 c write(iout,*) "f1",f1," f2",f2," hd1",hd1," hd2",hd2
431 c write (iout,*) "ssm",ssm," h1",h1," Ht",Ht," h2",h2
432 delta_inv=1.0d0/(xm-ssxm)
433 deltasq_inv=delta_inv*delta_inv
435 fac1=deltasq_inv*fac*(xm-rij)
436 fac2=deltasq_inv*fac*(rij-ssxm)
437 ed=delta_inv*(Ht*hd2-ssm*hd1)
439 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
440 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
441 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
442 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
445 ljm=-0.25D0*ljB*bb/aa
446 d_ljm(1)=-0.5D0*bb/aa*ljB
447 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
448 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
450 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
451 f1=(rij-ljxm)/(xm-ljxm)
452 f2=(rij-xm)/(ljxm-xm)
456 delta_inv=1.0d0/(ljxm-xm)
457 deltasq_inv=delta_inv*delta_inv
459 fac1=deltasq_inv*fac*(ljxm-rij)
460 fac2=deltasq_inv*fac*(rij-xm)
461 ed=delta_inv*(ljm*hd2-Ht*hd1)
463 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
464 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
465 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
466 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
468 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
470 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
476 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
477 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
478 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
480 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
481 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
482 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
483 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
486 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
488 c$$$ d_ljm(k)=ljm*d_ljB(k)
492 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
493 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
494 c$$$ d_ss(2)=akct*ssd
495 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
496 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
499 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
500 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
501 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
503 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
504 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
506 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
508 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
509 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
510 c$$$ h1=h_base(f1,hd1)
511 c$$$ h2=h_base(f2,hd2)
512 c$$$ eij=ss*h1+ljf*h2
513 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
514 c$$$ deltasq_inv=delta_inv*delta_inv
515 c$$$ fac=ljf*hd2-ss*hd1
516 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
517 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
518 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
519 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
520 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
521 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
522 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
524 c$$$ havebond=.false.
525 c$$$ if (ed.gt.0.0d0) havebond=.true.
526 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
533 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
534 c write(iout,'(a15,f12.2,f8.1,2i5)')
535 c & "SSBOND_E_FORM",totT,t_bath,i,j
539 dyn_ssbond_ij(i,j)=eij
540 else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
541 dyn_ssbond_ij(i,j)=1.0d300
544 c write(iout,'(a15,f12.2,f8.1,2i5)')
545 c & "SSBOND_E_BREAK",totT,t_bath,i,j
552 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
553 & "CHECKSTOP",rij,eij,ed
558 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
565 c-------END TESTING CODE
566 gg_lipi(3)=ssgradlipi*eij
567 gg_lipj(3)=ssgradlipj*eij
570 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
571 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
574 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
577 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
578 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
579 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
580 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
581 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
582 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
586 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
591 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
592 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
596 write (iout,*) "ssbond i",i," j",j," typ",itypi,itypj,
597 & " r",rij," ljxm",ljxm," ssxm",ssxm," xm",xm,
598 ^ " havebond",havebond," eij",eij
603 C-----------------------------------------------------------------------------
605 double precision function h_base(x,deriv)
606 c A smooth function going 0->1 in range [0,1]
607 c It should NOT be called outside range [0,1], it will not work there.
614 double precision deriv
620 c Two parabolas put together. First derivative zero at extrema
621 c$$$ if (x.lt.0.5D0) then
622 c$$$ h_base=2.0D0*x*x
626 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
627 c$$$ deriv=4.0D0*deriv
630 c Third degree polynomial. First derivative zero at extrema
631 h_base=x*x*(3.0d0-2.0d0*x)
632 deriv=6.0d0*x*(1.0d0-x)
634 c Fifth degree polynomial. First and second derivatives zero at extrema
636 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
638 c$$$ deriv=deriv*deriv
639 c$$$ deriv=30.0d0*xsq*deriv
644 c----------------------------------------------------------------------------
646 subroutine dyn_set_nss
647 c Adjust nss and other relevant variables based on dyn_ssbond_ij
655 include 'COMMON.SBRIDGE'
656 include 'COMMON.CHAIN'
657 include 'COMMON.IOUNITS'
658 c include 'COMMON.SETUP'
666 double precision emin
668 integer diff,allflag(maxdim),allnss,
669 & allihpb(maxdim),alljhpb(maxdim),
670 & newnss,newihpb(maxdim),newjhpb(maxdim)
672 c integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
673 c integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
678 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
687 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
691 if (allflag(i).eq.0 .and.
692 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
693 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
697 if (emin.lt.1.0d300) then
700 if (allflag(i).eq.0 .and.
701 & (allihpb(i).eq.allihpb(imin) .or.
702 & alljhpb(i).eq.allihpb(imin) .or.
703 & allihpb(i).eq.alljhpb(imin) .or.
704 & alljhpb(i).eq.alljhpb(imin))) then
711 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
715 if (allflag(i).eq.1) then
717 newihpb(newnss)=allihpb(i)
718 newjhpb(newnss)=alljhpb(i)
723 c if (nfgtasks.gt.1)then
725 c call MPI_Reduce(newnss,g_newnss,1,
726 c & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
727 c call MPI_Gather(newnss,1,MPI_INTEGER,
728 c & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
730 c do i=1,nfgtasks-1,1
731 c displ(i)=i_newnss(i-1)+displ(i-1)
733 c call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
734 c & g_newihpb,i_newnss,displ,MPI_INTEGER,
735 c & king,FG_COMM,IERR)
736 c call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
737 c & g_newjhpb,i_newnss,displ,MPI_INTEGER,
738 c & king,FG_COMM,IERR)
739 c if(fg_rank.eq.0) then
740 cc print *,'g_newnss',g_newnss
741 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
742 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
745 c newihpb(i)=g_newihpb(i)
746 c newjhpb(i)=g_newjhpb(i)
754 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
759 if (idssb(i).eq.newihpb(j) .and.
760 & jdssb(i).eq.newjhpb(j)) found=.true.
764 c if (.not.found.and.fg_rank.eq.0)
765 c & write(iout,'(a15,f12.2,f8.1,2i5)')
766 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
774 if (newihpb(i).eq.idssb(j) .and.
775 & newjhpb(i).eq.jdssb(j)) found=.true.
779 c if (.not.found.and.fg_rank.eq.0)
780 c & write(iout,'(a15,f12.2,f8.1,2i5)')
781 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
796 C-----------------------------------------------------------------------------
797 C-----------------------------------------------------------------------------
798 C-----------------------------------------------------------------------------
800 c$$$c-----------------------------------------------------------------------------
802 c$$$ subroutine ss_relax(i_in,j_in)
806 c$$$ include 'DIMENSIONS'
807 c$$$ include 'COMMON.VAR'
808 c$$$ include 'COMMON.CHAIN'
809 c$$$ include 'COMMON.IOUNITS'
810 c$$$ include 'COMMON.INTERACT'
812 c$$$c Input arguments
813 c$$$ integer i_in,j_in
815 c$$$c Local variables
816 c$$$ integer i,iretcode,nfun_sc
818 c$$$ double precision var(maxvar),e_sc,etot
825 c$$$ mask_side(i_in)=1
826 c$$$ mask_side(j_in)=1
828 c$$$c Minimize the two selected side-chains
829 c$$$ call overlap_sc(scfail) ! Better not fail!
830 c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
837 c$$$c-------------------------------------------------------------
839 c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
840 c$$$c Minimize side-chains only, starting from geom but without modifying
842 c$$$c If mask_r is already set, only the selected side-chains are minimized,
843 c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
847 c$$$ include 'DIMENSIONS'
848 c$$$ include 'COMMON.IOUNITS'
849 c$$$ include 'COMMON.VAR'
850 c$$$ include 'COMMON.CHAIN'
851 c$$$ include 'COMMON.GEO'
852 c$$$ include 'COMMON.MINIM'
854 c$$$ common /srutu/ icall
856 c$$$c Output arguments
857 c$$$ double precision etot_sc
858 c$$$ integer iretcode,nfun
860 c$$$c External functions/subroutines
861 c$$$ external func_sc,grad_sc,fdum
863 c$$$c Local variables
865 c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
867 c$$$ double precision rdum(1)
868 c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
870 c$$$ integer i,nvar_restr
873 c$$$cmc start_minim=.true.
874 c$$$ call deflt(2,iv,liv,lv,v)
875 c$$$* 12 means fresh start, dont call deflt
877 c$$$* max num of fun calls
878 c$$$ if (maxfun.eq.0) maxfun=500
880 c$$$* max num of iterations
881 c$$$ if (maxmin.eq.0) maxmin=1000
883 c$$$* controls output
885 c$$$* selects output unit
887 c$$$c iv(21)=iout ! DEBUG
888 c$$$c iv(21)=8 ! DEBUG
889 c$$$* 1 means to print out result
891 c$$$c iv(22)=1 ! DEBUG
892 c$$$* 1 means to print out summary stats
894 c$$$c iv(23)=1 ! DEBUG
895 c$$$* 1 means to print initial x and d
897 c$$$c iv(24)=1 ! DEBUG
898 c$$$* min val for v(radfac) default is 0.1
900 c$$$* max val for v(radfac) default is 4.0
903 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
904 c$$$* the sumsl default is 0.1
906 c$$$* false conv if (act fnctn decrease) .lt. v(34)
907 c$$$* the sumsl default is 100*machep
908 c$$$ v(34)=v(34)/100.0D0
909 c$$$* absolute convergence
910 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
912 c$$$* relative convergence
913 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
915 c$$$* controls initial step size
917 c$$$* large vals of d correspond to small components of step
921 c$$$ do i=nphi+1,nvar
925 c$$$ call geom_to_var(nvar,x)
926 c$$$ IF (mask_r) THEN
927 c$$$ do i=1,nres ! Just in case...
931 c$$$ call x2xx(x,xx,nvar_restr)
932 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
933 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
936 c$$$c When minimizing ALL side-chains, etotal_sc is a little
937 c$$$c faster if we don't set mask_r
943 c$$$ call x2xx(x,xx,nvar_restr)
944 c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
945 c$$$ & iv,liv,lv,v,idum,rdum,fdum)
948 c$$$ call var_to_geom(nvar,x)
949 c$$$ call chainbuild_sc
956 c$$$C--------------------------------------------------------------------------
958 c$$$ subroutine chainbuild_sc
960 c$$$ include 'DIMENSIONS'
961 c$$$ include 'COMMON.VAR'
962 c$$$ include 'COMMON.INTERACT'
964 c$$$c Local variables
969 c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
970 c$$$ call locate_side_chain(i)
977 c$$$C--------------------------------------------------------------------------
979 c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
983 c$$$ include 'DIMENSIONS'
984 c$$$ include 'COMMON.DERIV'
985 c$$$ include 'COMMON.VAR'
986 c$$$ include 'COMMON.MINIM'
987 c$$$ include 'COMMON.IOUNITS'
989 c$$$c Input arguments
991 c$$$ double precision x(maxvar)
992 c$$$ double precision ufparm
995 c$$$c Input/Output arguments
997 c$$$ integer uiparm(1)
998 c$$$ double precision urparm(1)
1000 c$$$c Output arguments
1001 c$$$ double precision f
1003 c$$$c Local variables
1004 c$$$ double precision energia(0:n_ene)
1006 c$$$c Variables used to intercept NaNs
1007 c$$$ double precision x_sum
1013 c$$$ icg=mod(nf,2)+1
1016 c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
1019 c$$$ x_sum=x_sum+x(i_NAN)
1021 c$$$c Calculate the energy only if the coordinates are ok
1022 c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
1023 c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
1029 c$$$ call var_to_geom_restr(n,x)
1031 c$$$ call chainbuild_sc
1032 c$$$ call etotal_sc(energia(0))
1034 c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
1043 c$$$c-------------------------------------------------------
1045 c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
1049 c$$$ include 'DIMENSIONS'
1050 c$$$ include 'COMMON.CHAIN'
1051 c$$$ include 'COMMON.DERIV'
1052 c$$$ include 'COMMON.VAR'
1053 c$$$ include 'COMMON.INTERACT'
1054 c$$$ include 'COMMON.MINIM'
1056 c$$$c Input arguments
1058 c$$$ double precision x(maxvar)
1059 c$$$ double precision ufparm
1060 c$$$ external ufparm
1062 c$$$c Input/Output arguments
1064 c$$$ integer uiparm(1)
1065 c$$$ double precision urparm(1)
1067 c$$$c Output arguments
1068 c$$$ double precision g(maxvar)
1070 c$$$c Local variables
1071 c$$$ double precision f,gphii,gthetai,galphai,gomegai
1072 c$$$ integer ig,ind,i,j,k,igall,ij
1075 c$$$ icg=mod(nf,2)+1
1076 c$$$ if (nf-nfl+1) 20,30,40
1077 c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1078 c$$$c write (iout,*) 'grad 20'
1079 c$$$ if (nf.eq.0) return
1081 c$$$ 30 call var_to_geom_restr(n,x)
1082 c$$$ call chainbuild_sc
1084 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1086 c$$$ 40 call cartder
1088 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1094 c$$$ IF (mask_phi(i+2).eq.1) THEN
1096 c$$$ do j=i+1,nres-1
1099 c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1100 c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1106 c$$$ ind=ind+nres-1-i
1113 c$$$ IF (mask_theta(i+2).eq.1) THEN
1116 c$$$ do j=i+1,nres-1
1119 c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1120 c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1125 c$$$ ind=ind+nres-1-i
1130 c$$$ if (itype(i).ne.10) then
1131 c$$$ IF (mask_side(i).eq.1) THEN
1135 c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1144 c$$$ if (itype(i).ne.10) then
1145 c$$$ IF (mask_side(i).eq.1) THEN
1149 c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1157 c$$$C Add the components corresponding to local energy terms.
1164 c$$$ if (mask_phi(i).eq.1) then
1166 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1172 c$$$ if (mask_theta(i).eq.1) then
1174 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1180 c$$$ if (itype(i).ne.10) then
1182 c$$$ if (mask_side(i).eq.1) then
1184 c$$$ g(ig)=g(ig)+gloc(igall,icg)
1191 c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1197 c$$$C-----------------------------------------------------------------------------
1199 c$$$ subroutine etotal_sc(energy_sc)
1203 c$$$ include 'DIMENSIONS'
1204 c$$$ include 'COMMON.VAR'
1205 c$$$ include 'COMMON.INTERACT'
1206 c$$$ include 'COMMON.DERIV'
1207 c$$$ include 'COMMON.FFIELD'
1209 c$$$c Output arguments
1210 c$$$ double precision energy_sc(0:n_ene)
1212 c$$$c Local variables
1213 c$$$ double precision evdw,escloc
1218 c$$$ energy_sc(i)=0.0D0
1221 c$$$ if (mask_r) then
1222 c$$$ call egb_sc(evdw)
1223 c$$$ call esc_sc(escloc)
1226 c$$$ call esc(escloc)
1229 c$$$ if (evdw.eq.1.0D20) then
1230 c$$$ energy_sc(0)=evdw
1232 c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
1234 c$$$ energy_sc(1)=evdw
1235 c$$$ energy_sc(12)=escloc
1238 c$$$C Sum up the components of the Cartesian gradient.
1242 c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
1249 c$$$C-----------------------------------------------------------------------------
1251 c$$$ subroutine egb_sc(evdw)
1253 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1254 c$$$C assuming the Gay-Berne potential of interaction.
1256 c$$$ implicit real*8 (a-h,o-z)
1257 c$$$ include 'DIMENSIONS'
1258 c$$$ include 'COMMON.GEO'
1259 c$$$ include 'COMMON.VAR'
1260 c$$$ include 'COMMON.LOCAL'
1261 c$$$ include 'COMMON.CHAIN'
1262 c$$$ include 'COMMON.DERIV'
1263 c$$$ include 'COMMON.NAMES'
1264 c$$$ include 'COMMON.INTERACT'
1265 c$$$ include 'COMMON.IOUNITS'
1266 c$$$ include 'COMMON.CALC'
1267 c$$$ include 'COMMON.CONTROL'
1270 c$$$ energy_dec=.false.
1271 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1274 c$$$c if (icall.eq.0) lprn=.false.
1276 c$$$ do i=iatsc_s,iatsc_e
1278 c$$$ itypi1=itype(i+1)
1282 c$$$ dxi=dc_norm(1,nres+i)
1283 c$$$ dyi=dc_norm(2,nres+i)
1284 c$$$ dzi=dc_norm(3,nres+i)
1285 c$$$c dsci_inv=dsc_inv(itypi)
1286 c$$$ dsci_inv=vbld_inv(i+nres)
1287 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1288 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1290 c$$$C Calculate SC interaction energy.
1292 c$$$ do iint=1,nint_gr(i)
1293 c$$$ do j=istart(i,iint),iend(i,iint)
1294 c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1297 c$$$c dscj_inv=dsc_inv(itypj)
1298 c$$$ dscj_inv=vbld_inv(j+nres)
1299 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1300 c$$$c & 1.0d0/vbld(j+nres)
1301 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1302 c$$$ sig0ij=sigma(itypi,itypj)
1303 c$$$ chi1=chi(itypi,itypj)
1304 c$$$ chi2=chi(itypj,itypi)
1305 c$$$ chi12=chi1*chi2
1306 c$$$ chip1=chip(itypi)
1307 c$$$ chip2=chip(itypj)
1308 c$$$ chip12=chip1*chip2
1309 c$$$ alf1=alp(itypi)
1310 c$$$ alf2=alp(itypj)
1311 c$$$ alf12=0.5D0*(alf1+alf2)
1312 c$$$C For diagnostics only!!!
1322 c$$$ xj=c(1,nres+j)-xi
1323 c$$$ yj=c(2,nres+j)-yi
1324 c$$$ zj=c(3,nres+j)-zi
1325 c$$$ dxj=dc_norm(1,nres+j)
1326 c$$$ dyj=dc_norm(2,nres+j)
1327 c$$$ dzj=dc_norm(3,nres+j)
1328 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1329 c$$$c write (iout,*) "j",j," dc_norm",
1330 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1331 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1332 c$$$ rij=dsqrt(rrij)
1333 c$$$C Calculate angle-dependent terms of energy and contributions to their
1335 c$$$ call sc_angular
1336 c$$$ sigsq=1.0D0/sigsq
1337 c$$$ sig=sig0ij*dsqrt(sigsq)
1338 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1339 c$$$c for diagnostics; uncomment
1340 c$$$c rij_shift=1.2*sig0ij
1341 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1342 c$$$ if (rij_shift.le.0.0D0) then
1344 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1345 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1346 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1349 c$$$ sigder=-sig*sigsq
1350 c$$$c---------------------------------------------------------------
1351 c$$$ rij_shift=1.0D0/rij_shift
1352 c$$$ fac=rij_shift**expon
1353 c$$$ e1=fac*fac*aa(itypi,itypj)
1354 c$$$ e2=fac*bb(itypi,itypj)
1355 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1356 c$$$ eps2der=evdwij*eps3rt
1357 c$$$ eps3der=evdwij*eps2rt
1358 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1359 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1360 c$$$ evdwij=evdwij*eps2rt*eps3rt
1361 c$$$ evdw=evdw+evdwij
1363 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1364 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1365 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1366 c$$$ & restyp(itypi),i,restyp(itypj),j,
1367 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1368 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1369 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1373 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1374 c$$$ & 'evdw',i,j,evdwij
1376 c$$$C Calculate gradient components.
1377 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1378 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1379 c$$$ sigder=fac*sigder
1382 c$$$C Calculate the radial part of the gradient
1386 c$$$C Calculate angular part of the gradient.
1392 c$$$ energy_dec=.false.
1396 c$$$c-----------------------------------------------------------------------------
1398 c$$$ subroutine esc_sc(escloc)
1399 c$$$C Calculate the local energy of a side chain and its derivatives in the
1400 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
1401 c$$$C ALPHA and OMEGA.
1402 c$$$ implicit real*8 (a-h,o-z)
1403 c$$$ include 'DIMENSIONS'
1404 c$$$ include 'COMMON.GEO'
1405 c$$$ include 'COMMON.LOCAL'
1406 c$$$ include 'COMMON.VAR'
1407 c$$$ include 'COMMON.INTERACT'
1408 c$$$ include 'COMMON.DERIV'
1409 c$$$ include 'COMMON.CHAIN'
1410 c$$$ include 'COMMON.IOUNITS'
1411 c$$$ include 'COMMON.NAMES'
1412 c$$$ include 'COMMON.FFIELD'
1413 c$$$ include 'COMMON.CONTROL'
1414 c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1415 c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
1416 c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
1417 c$$$ delta=0.02d0*pi
1419 c$$$c write (iout,'(a)') 'ESC'
1420 c$$$ do i=loc_start,loc_end
1421 c$$$ IF (mask_side(i).eq.1) THEN
1423 c$$$ if (it.eq.10) goto 1
1424 c$$$ nlobit=nlob(it)
1425 c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
1426 c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1427 c$$$ theti=theta(i+1)-pipol
1428 c$$$ x(1)=dtan(theti)
1432 c$$$ if (x(2).gt.pi-delta) then
1434 c$$$ xtemp(2)=pi-delta
1436 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1438 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1439 c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1440 c$$$ & escloci,dersc(2))
1441 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1442 c$$$ & ddersc0(1),dersc(1))
1443 c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1444 c$$$ & ddersc0(3),dersc(3))
1445 c$$$ xtemp(2)=pi-delta
1446 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1448 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1449 c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1450 c$$$ & dersc0(2),esclocbi,dersc02)
1451 c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1452 c$$$ & dersc12,dersc01)
1453 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1454 c$$$ dersc0(1)=dersc01
1455 c$$$ dersc0(2)=dersc02
1456 c$$$ dersc0(3)=0.0d0
1458 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1460 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1461 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1462 c$$$c & esclocbi,ss,ssd
1463 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1464 c$$$c escloci=esclocbi
1465 c$$$c write (iout,*) escloci
1466 c$$$ else if (x(2).lt.delta) then
1470 c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1472 c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1473 c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1474 c$$$ & escloci,dersc(2))
1475 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1476 c$$$ & ddersc0(1),dersc(1))
1477 c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1478 c$$$ & ddersc0(3),dersc(3))
1480 c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1482 c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1483 c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1484 c$$$ & dersc0(2),esclocbi,dersc02)
1485 c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1486 c$$$ & dersc12,dersc01)
1487 c$$$ dersc0(1)=dersc01
1488 c$$$ dersc0(2)=dersc02
1489 c$$$ dersc0(3)=0.0d0
1490 c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
1492 c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1494 c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1495 c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1496 c$$$c & esclocbi,ss,ssd
1497 c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
1498 c$$$c write (iout,*) escloci
1500 c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
1503 c$$$ escloc=escloc+escloci
1504 c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1505 c$$$ & 'escloc',i,escloci
1506 c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1508 c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1509 c$$$ & wscloc*dersc(1)
1510 c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
1511 c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1518 c$$$C-----------------------------------------------------------------------------
1520 c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
1522 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1523 c$$$C assuming the Gay-Berne potential of interaction.
1525 c$$$ implicit real*8 (a-h,o-z)
1526 c$$$ include 'DIMENSIONS'
1527 c$$$ include 'COMMON.GEO'
1528 c$$$ include 'COMMON.VAR'
1529 c$$$ include 'COMMON.LOCAL'
1530 c$$$ include 'COMMON.CHAIN'
1531 c$$$ include 'COMMON.DERIV'
1532 c$$$ include 'COMMON.NAMES'
1533 c$$$ include 'COMMON.INTERACT'
1534 c$$$ include 'COMMON.IOUNITS'
1535 c$$$ include 'COMMON.CALC'
1536 c$$$ include 'COMMON.CONTROL'
1539 c$$$ energy_dec=.false.
1540 c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1544 c$$$c$$$ do i=iatsc_s,iatsc_e
1547 c$$$ itypi1=itype(i+1)
1551 c$$$ dxi=dc_norm(1,nres+i)
1552 c$$$ dyi=dc_norm(2,nres+i)
1553 c$$$ dzi=dc_norm(3,nres+i)
1554 c$$$c dsci_inv=dsc_inv(itypi)
1555 c$$$ dsci_inv=vbld_inv(i+nres)
1556 c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1557 c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1559 c$$$C Calculate SC interaction energy.
1561 c$$$c$$$ do iint=1,nint_gr(i)
1562 c$$$c$$$ do j=istart(i,iint),iend(i,iint)
1566 c$$$c dscj_inv=dsc_inv(itypj)
1567 c$$$ dscj_inv=vbld_inv(j+nres)
1568 c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1569 c$$$c & 1.0d0/vbld(j+nres)
1570 c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1571 c$$$ sig0ij=sigma(itypi,itypj)
1572 c$$$ chi1=chi(itypi,itypj)
1573 c$$$ chi2=chi(itypj,itypi)
1574 c$$$ chi12=chi1*chi2
1575 c$$$ chip1=chip(itypi)
1576 c$$$ chip2=chip(itypj)
1577 c$$$ chip12=chip1*chip2
1578 c$$$ alf1=alp(itypi)
1579 c$$$ alf2=alp(itypj)
1580 c$$$ alf12=0.5D0*(alf1+alf2)
1581 c$$$C For diagnostics only!!!
1591 c$$$ xj=c(1,nres+j)-xi
1592 c$$$ yj=c(2,nres+j)-yi
1593 c$$$ zj=c(3,nres+j)-zi
1594 c$$$ dxj=dc_norm(1,nres+j)
1595 c$$$ dyj=dc_norm(2,nres+j)
1596 c$$$ dzj=dc_norm(3,nres+j)
1597 c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1598 c$$$c write (iout,*) "j",j," dc_norm",
1599 c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1600 c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1601 c$$$ rij=dsqrt(rrij)
1602 c$$$C Calculate angle-dependent terms of energy and contributions to their
1604 c$$$ call sc_angular
1605 c$$$ sigsq=1.0D0/sigsq
1606 c$$$ sig=sig0ij*dsqrt(sigsq)
1607 c$$$ rij_shift=1.0D0/rij-sig+sig0ij
1608 c$$$c for diagnostics; uncomment
1609 c$$$c rij_shift=1.2*sig0ij
1610 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1611 c$$$ if (rij_shift.le.0.0D0) then
1613 c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1614 c$$$cd & restyp(itypi),i,restyp(itypj),j,
1615 c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1618 c$$$ sigder=-sig*sigsq
1619 c$$$c---------------------------------------------------------------
1620 c$$$ rij_shift=1.0D0/rij_shift
1621 c$$$ fac=rij_shift**expon
1622 c$$$ e1=fac*fac*aa(itypi,itypj)
1623 c$$$ e2=fac*bb(itypi,itypj)
1624 c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1625 c$$$ eps2der=evdwij*eps3rt
1626 c$$$ eps3der=evdwij*eps2rt
1627 c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1628 c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1629 c$$$ evdwij=evdwij*eps2rt*eps3rt
1630 c$$$ evdw=evdw+evdwij
1632 c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1633 c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1634 c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1635 c$$$ & restyp(itypi),i,restyp(itypj),j,
1636 c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
1637 c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1638 c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1642 c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
1643 c$$$ & 'evdw',i,j,evdwij
1645 c$$$C Calculate gradient components.
1646 c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
1647 c$$$ fac=-expon*(e1+evdwij)*rij_shift
1648 c$$$ sigder=fac*sigder
1651 c$$$C Calculate the radial part of the gradient
1655 c$$$C Calculate angular part of the gradient.
1658 c$$$c$$$ enddo ! iint
1660 c$$$ energy_dec=.false.
1664 c$$$C-----------------------------------------------------------------------------
1666 c$$$ subroutine perturb_side_chain(i,angle)
1670 c$$$ include 'DIMENSIONS'
1671 c$$$ include 'COMMON.CHAIN'
1672 c$$$ include 'COMMON.GEO'
1673 c$$$ include 'COMMON.VAR'
1674 c$$$ include 'COMMON.LOCAL'
1675 c$$$ include 'COMMON.IOUNITS'
1677 c$$$c External functions
1678 c$$$ external ran_number
1679 c$$$ double precision ran_number
1681 c$$$c Input arguments
1683 c$$$ double precision angle ! In degrees
1685 c$$$c Local variables
1687 c$$$ double precision rad_ang,rand_v(3),length,cost,sint
1691 c$$$ rad_ang=angle*deg2rad
1694 c$$$ do while (length.lt.0.01)
1695 c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
1696 c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
1697 c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
1698 c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1699 c$$$ + rand_v(3)*rand_v(3)
1700 c$$$ length=sqrt(length)
1701 c$$$ rand_v(1)=rand_v(1)/length
1702 c$$$ rand_v(2)=rand_v(2)/length
1703 c$$$ rand_v(3)=rand_v(3)/length
1704 c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1705 c$$$ + rand_v(3)*dc_norm(3,i_sc)
1706 c$$$ length=1.0D0-cost*cost
1707 c$$$ if (length.lt.0.0D0) length=0.0D0
1708 c$$$ length=sqrt(length)
1709 c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1710 c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1711 c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1713 c$$$ rand_v(1)=rand_v(1)/length
1714 c$$$ rand_v(2)=rand_v(2)/length
1715 c$$$ rand_v(3)=rand_v(3)/length
1717 c$$$ cost=dcos(rad_ang)
1718 c$$$ sint=dsin(rad_ang)
1719 c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1720 c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1721 c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1722 c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1723 c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1724 c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1725 c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
1726 c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
1727 c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
1729 c$$$ call chainbuild_cart
1734 c$$$c----------------------------------------------------------------------------
1736 c$$$ subroutine ss_relax3(i_in,j_in)
1740 c$$$ include 'DIMENSIONS'
1741 c$$$ include 'COMMON.VAR'
1742 c$$$ include 'COMMON.CHAIN'
1743 c$$$ include 'COMMON.IOUNITS'
1744 c$$$ include 'COMMON.INTERACT'
1746 c$$$c External functions
1747 c$$$ external ran_number
1748 c$$$ double precision ran_number
1750 c$$$c Input arguments
1751 c$$$ integer i_in,j_in
1753 c$$$c Local variables
1754 c$$$ double precision energy_sc(0:n_ene),etot
1755 c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
1756 c$$$ double precision ang_pert,rand_fact,exp_fact,beta
1757 c$$$ integer n,i_pert,i
1758 c$$$ logical notdone
1767 c$$$ mask_side(i_in)=1
1768 c$$$ mask_side(j_in)=1
1770 c$$$ call etotal_sc(energy_sc)
1771 c$$$ etot=energy_sc(0)
1772 c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
1773 c$$$c + energy_sc(1),energy_sc(12)
1777 c$$$ do while (notdone)
1778 c$$$ if (mod(n,2).eq.0) then
1786 c$$$ org_dc(i)=dc(i,i_pert+nres)
1787 c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
1788 c$$$ org_c(i)=c(i,i_pert+nres)
1790 c$$$ ang_pert=ran_number(0.0D0,3.0D0)
1791 c$$$ call perturb_side_chain(i_pert,ang_pert)
1792 c$$$ call etotal_sc(energy_sc)
1793 c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
1794 c$$$ rand_fact=ran_number(0.0D0,1.0D0)
1795 c$$$ if (rand_fact.lt.exp_fact) then
1796 c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
1797 c$$$c + energy_sc(1),energy_sc(12)
1798 c$$$ etot=energy_sc(0)
1800 c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
1801 c$$$c + energy_sc(1),energy_sc(12)
1803 c$$$ dc(i,i_pert+nres)=org_dc(i)
1804 c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
1805 c$$$ c(i,i_pert+nres)=org_c(i)
1809 c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1817 c$$$c----------------------------------------------------------------------------
1819 c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1821 c$$$ include 'DIMENSIONS'
1823 c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
1824 c$$$*********************************************************************
1825 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1826 c$$$* the calling subprogram. *
1827 c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
1828 c$$$* calculated in the usual pythagorean way. *
1829 c$$$* absolute convergence occurs when the function is within v(31) of *
1830 c$$$* zero. unless you know the minimum value in advance, abs convg *
1831 c$$$* is probably not useful. *
1832 c$$$* relative convergence is when the model predicts that the function *
1833 c$$$* will decrease by less than v(32)*abs(fun). *
1834 c$$$*********************************************************************
1835 c$$$ include 'COMMON.IOUNITS'
1836 c$$$ include 'COMMON.VAR'
1837 c$$$ include 'COMMON.GEO'
1838 c$$$ include 'COMMON.MINIM'
1839 c$$$ include 'COMMON.CHAIN'
1841 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1842 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1843 c$$$ + orig_ss_dist(maxres2,maxres2)
1845 c$$$ double precision etot
1846 c$$$ integer iretcode,nfun,i_in,j_in
1849 c$$$ double precision dist
1850 c$$$ external ss_func,fdum
1851 c$$$ double precision ss_func,fdum
1853 c$$$ integer iv(liv),uiparm(2)
1854 c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
1858 c$$$ call deflt(2,iv,liv,lv,v)
1859 c$$$* 12 means fresh start, dont call deflt
1861 c$$$* max num of fun calls
1862 c$$$ if (maxfun.eq.0) maxfun=500
1864 c$$$* max num of iterations
1865 c$$$ if (maxmin.eq.0) maxmin=1000
1867 c$$$* controls output
1869 c$$$* selects output unit
1872 c$$$* 1 means to print out result
1874 c$$$* 1 means to print out summary stats
1876 c$$$* 1 means to print initial x and d
1878 c$$$* min val for v(radfac) default is 0.1
1880 c$$$* max val for v(radfac) default is 4.0
1883 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1884 c$$$* the sumsl default is 0.1
1886 c$$$* false conv if (act fnctn decrease) .lt. v(34)
1887 c$$$* the sumsl default is 100*machep
1888 c$$$ v(34)=v(34)/100.0D0
1889 c$$$* absolute convergence
1890 c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
1893 c$$$* relative convergence
1894 c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
1897 c$$$* controls initial step size
1899 c$$$* large vals of d correspond to small components of step
1906 c$$$ orig_ss_dc(j,i)=dc(j,i)
1909 c$$$ call geom_to_var(nvar,orig_ss_var)
1913 c$$$ orig_ss_dist(j,i)=dist(j,i)
1914 c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
1915 c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
1916 c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1928 c$$$ if (ialph(i,1).gt.0) then
1931 c$$$ x(k)=dc(j,i+nres)
1938 c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1941 c$$$ nfun=iv(6)+iv(30)
1951 c$$$ if (ialph(i,1).gt.0) then
1954 c$$$ dc(j,i+nres)=x(k)
1958 c$$$ call chainbuild_cart
1963 c$$$C-----------------------------------------------------------------------------
1965 c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
1967 c$$$ include 'DIMENSIONS'
1968 c$$$ include 'COMMON.DERIV'
1969 c$$$ include 'COMMON.IOUNITS'
1970 c$$$ include 'COMMON.VAR'
1971 c$$$ include 'COMMON.CHAIN'
1972 c$$$ include 'COMMON.INTERACT'
1973 c$$$ include 'COMMON.SBRIDGE'
1975 c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1976 c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1977 c$$$ + orig_ss_dist(maxres2,maxres2)
1980 c$$$ double precision x(maxres6)
1982 c$$$ double precision f
1983 c$$$ integer uiparm(2)
1984 c$$$ real*8 urparm(1)
1985 c$$$ external ufparm
1986 c$$$ double precision ufparm
1989 c$$$ double precision dist
1991 c$$$ integer i,j,k,ss_i,ss_j
1992 c$$$ double precision tempf,var(maxvar)
2007 c$$$ if (ialph(i,1).gt.0) then
2010 c$$$ dc(j,i+nres)=x(k)
2014 c$$$ call chainbuild_cart
2016 c$$$ call geom_to_var(nvar,var)
2018 c$$$c Constraints on all angles
2020 c$$$ tempf=var(i)-orig_ss_var(i)
2021 c$$$ f=f+tempf*tempf
2024 c$$$c Constraints on all distances
2026 c$$$ if (i.gt.1) then
2027 c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
2028 c$$$ f=f+tempf*tempf
2031 c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
2032 c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
2033 c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
2034 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2035 c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
2036 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2037 c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
2038 c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
2042 c$$$c Constraints for the relevant CYS-CYS
2043 c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
2044 c$$$ f=f+tempf*tempf
2045 c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
2047 c$$$c$$$ if (nf.ne.nfl) then
2048 c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
2049 c$$$c$$$ + f,dist(5+nres,14+nres)
2057 c$$$C-----------------------------------------------------------------------------
2058 c$$$C-----------------------------------------------------------------------------
2059 subroutine triple_ssbond_ene(resi,resj,resk,eij)
2060 include 'DIMENSIONS'
2061 include 'COMMON.SBRIDGE'
2062 include 'COMMON.CHAIN'
2063 include 'COMMON.DERIV'
2064 include 'COMMON.LOCAL'
2065 include 'COMMON.INTERACT'
2066 include 'COMMON.VAR'
2067 include 'COMMON.IOUNITS'
2068 include 'COMMON.CALC'
2075 c External functions
2076 double precision h_base
2080 integer resi,resj,resk
2083 double precision eij,eij1,eij2,eij3
2087 c integer itypi,itypj,k,l
2088 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2089 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2090 double precision xik,yik,zik,xjk,yjk,zjk
2091 double precision sig0ij,ljd,sig,fac,e1,e2
2092 double precision dcosom1(3),dcosom2(3),ed
2093 double precision pom1,pom2
2094 double precision ljA,ljB,ljXs
2095 double precision d_ljB(1:3)
2096 double precision ssA,ssB,ssC,ssXs
2097 double precision ssxm,ljxm,ssm,ljm
2098 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2099 if (dtriss.eq.0) return
2103 C write(iout,*) resi,resj,resk
2105 dxi=dc_norm(1,nres+i)
2106 dyi=dc_norm(2,nres+i)
2107 dzi=dc_norm(3,nres+i)
2108 dsci_inv=vbld_inv(i+nres)
2118 dxj=dc_norm(1,nres+j)
2119 dyj=dc_norm(2,nres+j)
2120 dzj=dc_norm(3,nres+j)
2121 dscj_inv=vbld_inv(j+nres)
2127 dxk=dc_norm(1,nres+k)
2128 dyk=dc_norm(2,nres+k)
2129 dzk=dc_norm(3,nres+k)
2130 dscj_inv=vbld_inv(k+nres)
2140 rrij=(xij*xij+yij*yij+zij*zij)
2141 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
2142 rrik=(xik*xik+yik*yik+zik*zik)
2144 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2146 C there are three combination of distances for each trisulfide bonds
2147 C The first case the ith atom is the center
2148 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2149 C distance y is second distance the a,b,c,d are parameters derived for
2150 C this problem d parameter was set as a penalty currenlty set to 1.
2151 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2152 C second case jth atom is center
2153 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2154 C the third case kth atom is the center
2155 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2160 C write(iout,*)i,j,k,eij
2161 C The energy penalty calculated now time for the gradient part
2162 C derivative over rij
2163 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2164 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
2169 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2170 gvdwx(m,j)=gvdwx(m,j)+gg(m)
2173 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2174 gvdwc(l,j)=gvdwc(l,j)+gg(l)
2176 C now derivative over rik
2177 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2178 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2183 gvdwx(m,i)=gvdwx(m,i)-gg(m)
2184 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2187 gvdwc(l,i)=gvdwc(l,i)-gg(l)
2188 gvdwc(l,k)=gvdwc(l,k)+gg(l)
2190 C now derivative over rjk
2191 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2192 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2197 gvdwx(m,j)=gvdwx(m,j)-gg(m)
2198 gvdwx(m,k)=gvdwx(m,k)+gg(m)
2201 gvdwc(l,j)=gvdwc(l,j)-gg(l)
2202 gvdwc(l,k)=gvdwc(l,k)+gg(l)