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-----------------------------------------------------------------------------
85 subroutine dyn_ssbond_ene(resi,resj,eij)
88 include 'DIMENSIONS.ZSCOPT'
89 include 'COMMON.SBRIDGE'
90 include 'COMMON.CHAIN'
91 include 'COMMON.DERIV'
92 include 'COMMON.LOCAL'
93 include 'COMMON.INTERACT'
95 include 'COMMON.IOUNITS'
97 include 'COMMON.NAMES'
105 double precision h_base
116 c integer itypi,itypj,k,l
117 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
118 double precision sig0ij,ljd,sig,fac,e1,e2
119 double precision dcosom1(3),dcosom2(3),ed
120 double precision pom1,pom2
121 double precision ljA,ljB,ljXs
122 double precision d_ljB(1:3)
123 double precision ssA,ssB,ssC,ssXs
124 double precision ssxm,ljxm,ssm,ljm
125 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
126 double precision f1,f2,h1,h2,hd1,hd2
127 double precision omega,delta_inv,deltasq_inv,fac1,fac2
129 double precision xm,d_xm(1:3)
130 double precision sslipi,sslipj,ssgradlipi,ssgradlipj
131 integer ici,icj,itypi,itypj
132 double precision boxshift,sscale,sscagrad
133 double precision aa,bb
134 c-------END FIRST METHOD
135 c-------SECOND METHOD
136 c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
137 c-------END SECOND METHOD
140 logical checkstop,transgrad
141 common /sschecks/ checkstop,transgrad
143 integer icheck,nicheck,jcheck,njcheck
144 double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
145 c-------END TESTING CODE
152 if (ici.eq.0 .or. icj.eq.0) then
153 write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
154 & "Attempt to create",
155 & " a disulfide link between non-cysteine residues ",restyp(i),i,
160 dxi=dc_norm(1,nres+i)
161 dyi=dc_norm(2,nres+i)
162 dzi=dc_norm(3,nres+i)
163 dsci_inv=vbld_inv(i+nres)
167 call to_box(xi,yi,zi)
168 C define scaling factor for lipids
170 C if (positi.le.0) positi=positi+boxzsize
172 C first for peptide groups
173 c for each residue check if it is in lipid or lipid water border area
174 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
179 call to_box(xj,yj,zj)
180 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
181 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
182 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
183 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
184 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
185 xj=boxshift(xj-xi,boxxsize)
186 yj=boxshift(yj-yi,boxysize)
187 zj=boxshift(zj-zi,boxzsize)
188 dxj=dc_norm(1,nres+j)
189 dyj=dc_norm(2,nres+j)
190 dzj=dc_norm(3,nres+j)
191 dscj_inv=vbld_inv(j+nres)
193 chi1=chi(itypi,itypj)
194 chi2=chi(itypj,itypi)
201 alf12=0.5D0*(alf1+alf2)
203 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
204 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
205 sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
206 sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
207 c The following are set in sc_angular
211 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
212 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
213 c om12=dxi*dxj+dyi*dyj+dzi*dzj
215 rij=1.0D0/rij ! Reset this so it makes sense
217 sig0ij=sigma(itypi,itypj)
218 sig=sig0ij*dsqrt(1.0D0/sigsq)
221 ljA=eps1*eps2rt**2*eps3rt**2
224 ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
229 deltat12=om2-om1+2.0d0
234 & +akth*(deltat1*deltat1+deltat2*deltat2)
235 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
236 ssxm=ssXs-0.5D0*ssB/ssA
239 c$$$c Some extra output
240 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
241 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
242 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
243 c$$$ if (ssx0.gt.0.0d0) then
244 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
248 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
249 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
250 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
252 c-------END TESTING CODE
255 c Stop and plot energy and derivative as a function of distance
257 ssm=ssC-0.25D0*ssB*ssB/ssA
258 ljm=-0.25D0*ljB*bb/aa
260 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
268 if (.not.checkstop) then
275 if (checkstop) rij=(ssxm-1.0d0)+
276 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
277 c-------END TESTING CODE
279 if (rij.gt.ljxm) then
282 fac=(1.0D0/ljd)**expon
285 eij=eps1*eps2rt*eps3rt*(e1+e2)
288 eij=eij*eps2rt*eps3rt*sss
291 e1=e1*eps1*eps2rt**2*eps3rt**2
292 ed=-expon*(e1+eij)/ljd
294 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
295 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
296 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
297 eom12=eij*eps1_om12+eps2der*eps2rt_om12
298 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
299 else if (rij.lt.ssxm) then
302 eij=ssA*ssd*ssd+ssB*ssd+ssC
304 ed=2*akcm*ssd+akct*deltat12
305 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
307 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
308 eom1=-2*akth*deltat1-pom1-om2*pom2
309 eom2= 2*akth*deltat2+pom1-om1*pom2
312 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
314 d_ssxm(1)=0.5D0*akct/ssA
318 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
319 d_ljxm(2)=d_ljxm(1)*sigsq_om2
320 d_ljxm(3)=d_ljxm(1)*sigsq_om12
321 d_ljxm(1)=d_ljxm(1)*sigsq_om1
323 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
326 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
330 ssm=ssC-0.25D0*ssB*ssB/ssA
331 d_ssm(1)=0.5D0*akct*ssB/ssA
332 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
333 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
335 f1=(rij-xm)/(ssxm-xm)
336 f2=(rij-ssxm)/(xm-ssxm)
340 delta_inv=1.0d0/(xm-ssxm)
341 deltasq_inv=delta_inv*delta_inv
343 fac1=deltasq_inv*fac*(xm-rij)
344 fac2=deltasq_inv*fac*(rij-ssxm)
345 ed=delta_inv*(Ht*hd2-ssm*hd1)
347 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
348 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
349 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
350 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
353 ljm=-0.25D0*ljB*bb/aa
354 d_ljm(1)=-0.5D0*bb/aa*ljB
355 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
356 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
358 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
359 f1=(rij-ljxm)/(xm-ljxm)
360 f2=(rij-xm)/(ljxm-xm)
364 delta_inv=1.0d0/(ljxm-xm)
365 deltasq_inv=delta_inv*delta_inv
367 fac1=deltasq_inv*fac*(ljxm-rij)
368 fac2=deltasq_inv*fac*(rij-xm)
369 ed=delta_inv*(ljm*hd2-Ht*hd1)
371 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
372 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
373 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
374 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
376 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
378 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
384 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
385 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
386 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
388 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
389 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
390 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
391 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
394 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
396 c$$$ d_ljm(k)=ljm*d_ljB(k)
400 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
401 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
402 c$$$ d_ss(2)=akct*ssd
403 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
404 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
407 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
408 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
409 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
411 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
412 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
414 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
416 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
417 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
418 c$$$ h1=h_base(f1,hd1)
419 c$$$ h2=h_base(f2,hd2)
420 c$$$ eij=ss*h1+ljf*h2
421 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
422 c$$$ deltasq_inv=delta_inv*delta_inv
423 c$$$ fac=ljf*hd2-ss*hd1
424 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
425 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
426 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
427 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
428 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
429 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
430 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
432 c$$$ havebond=.false.
433 c$$$ if (ed.gt.0.0d0) havebond=.true.
434 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
441 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
442 c write(iout,'(a15,f12.2,f8.1,2i5)')
443 c & "SSBOND_E_FORM",totT,t_bath,i,j
447 dyn_ssbond_ij(ici,icj)=eij
448 else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300)
450 dyn_ssbond_ij(ici,icj)=1.0d300
453 c write(iout,'(a15,f12.2,f8.1,2i5)')
454 c & "SSBOND_E_BREAK",totT,t_bath,i,j
461 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
462 & "CHECKSTOP",rij,eij,ed
467 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
474 c-------END TESTING CODE
475 gg_lipi(3)=ssgradlipi*eij
476 gg_lipj(3)=ssgradlipj*eij
479 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
480 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
483 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
486 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
487 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
488 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
489 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
490 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
491 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
495 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
500 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
501 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
506 C-----------------------------------------------------------------------------
508 double precision function h_base(x,deriv)
509 c A smooth function going 0->1 in range [0,1]
510 c It should NOT be called outside range [0,1], it will not work there.
517 double precision deriv
523 c Two parabolas put together. First derivative zero at extrema
524 c$$$ if (x.lt.0.5D0) then
525 c$$$ h_base=2.0D0*x*x
529 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
530 c$$$ deriv=4.0D0*deriv
533 c Third degree polynomial. First derivative zero at extrema
534 h_base=x*x*(3.0d0-2.0d0*x)
535 deriv=6.0d0*x*(1.0d0-x)
537 c Fifth degree polynomial. First and second derivatives zero at extrema
539 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
541 c$$$ deriv=deriv*deriv
542 c$$$ deriv=30.0d0*xsq*deriv
546 c----------------------------------------------------------------------------
547 subroutine dyn_set_nss
548 c Adjust nss and other relevant variables based on dyn_ssbond_ij
556 include 'COMMON.SBRIDGE'
557 include 'COMMON.CHAIN'
558 include 'COMMON.IOUNITS'
559 C include 'COMMON.SETUP'
562 C include 'COMMON.MD'
567 double precision emin
569 integer diff,allflag(maxdim_cont),allnss,
570 & allihpb(maxdim_cont),alljhpb(maxdim_cont),
571 & newnss,newihpb(maxdim_cont),newjhpb(maxdim_cont)
573 integer i_newnss(1024),displ(0:1024)
574 integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss
579 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
588 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
592 if (allflag(i).eq.0 .and.
593 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
594 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
598 if (emin.lt.1.0d300) then
601 if (allflag(i).eq.0 .and.
602 & (allihpb(i).eq.allihpb(imin) .or.
603 & alljhpb(i).eq.allihpb(imin) .or.
604 & allihpb(i).eq.alljhpb(imin) .or.
605 & alljhpb(i).eq.alljhpb(imin))) then
612 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
616 if (allflag(i).eq.1) then
618 newihpb(newnss)=allihpb(i)
619 newjhpb(newnss)=alljhpb(i)
624 if (nfgtasks.gt.1)then
626 call MPI_Reduce(newnss,g_newnss,1,
627 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
628 call MPI_Gather(newnss,1,MPI_INTEGER,
629 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
632 displ(i)=i_newnss(i-1)+displ(i-1)
634 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
635 & g_newihpb,i_newnss,displ,MPI_INTEGER,
637 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
638 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
640 if(fg_rank.eq.0) then
641 c print *,'g_newnss',g_newnss
642 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
643 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
646 newihpb(i)=g_newihpb(i)
647 newjhpb(i)=g_newjhpb(i)
655 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
660 if (idssb(i).eq.newihpb(j) .and.
661 & jdssb(i).eq.newjhpb(j)) found=.true.
665 c if (.not.found.and.fg_rank.eq.0)
666 c & write(iout,'(a15,f12.2,f8.1,2i5)')
667 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
675 if (newihpb(i).eq.idssb(j) .and.
676 & newjhpb(i).eq.jdssb(j)) found=.true.
680 c if (.not.found.and.fg_rank.eq.0)
681 c & write(iout,'(a15,f12.2,f8.1,2i5)')
682 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
696 c----------------------------------------------------------------------------
700 subroutine read_ssHist
705 include "DIMENSIONS.FREE"
706 include 'COMMON.FREE'
710 character*80 controlcard
713 call card_concat(controlcard,.true.)
715 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
722 c$$$C----------------------------------------------------------------------------
723 subroutine triple_ssbond_ene(resi,resj,resk,eij)
725 include 'COMMON.SBRIDGE'
726 include 'COMMON.CHAIN'
727 include 'COMMON.DERIV'
728 include 'COMMON.LOCAL'
729 include 'COMMON.INTERACT'
731 include 'COMMON.IOUNITS'
732 include 'COMMON.CALC'
735 C include 'COMMON.MD'
740 double precision h_base
744 integer resi,resj,resk
747 double precision eij,eij1,eij2,eij3
751 c integer itypi,itypj,k,l
752 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
753 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
754 double precision xik,yik,zik,xjk,yjk,zjk
755 double precision sig0ij,ljd,sig,fac,e1,e2
756 double precision dcosom1(3),dcosom2(3),ed
757 double precision pom1,pom2
758 double precision ljA,ljB,ljXs
759 double precision d_ljB(1:3)
760 double precision ssA,ssB,ssC,ssXs
761 double precision ssxm,ljxm,ssm,ljm
762 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
767 C write(iout,*) resi,resj,resk
769 dxi=dc_norm(1,nres+i)
770 dyi=dc_norm(2,nres+i)
771 dzi=dc_norm(3,nres+i)
772 dsci_inv=vbld_inv(i+nres)
782 dxj=dc_norm(1,nres+j)
783 dyj=dc_norm(2,nres+j)
784 dzj=dc_norm(3,nres+j)
785 dscj_inv=vbld_inv(j+nres)
791 dxk=dc_norm(1,nres+k)
792 dyk=dc_norm(2,nres+k)
793 dzk=dc_norm(3,nres+k)
794 dscj_inv=vbld_inv(k+nres)
804 rrij=(xij*xij+yij*yij+zij*zij)
805 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
806 rrik=(xik*xik+yik*yik+zik*zik)
808 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
810 C there are three combination of distances for each trisulfide bonds
811 C The first case the ith atom is the center
812 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
813 C distance y is second distance the a,b,c,d are parameters derived for
814 C this problem d parameter was set as a penalty currenlty set to 1.
815 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
816 C second case jth atom is center
817 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
818 C the third case kth atom is the center
819 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
824 C write(iout,*)i,j,k,eij
825 C The energy penalty calculated now time for the gradient part
826 C derivative over rij
827 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
828 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
833 gvdwx(m,i)=gvdwx(m,i)-gg(m)
834 gvdwx(m,j)=gvdwx(m,j)+gg(m)
837 gvdwc(l,i)=gvdwc(l,i)-gg(l)
838 gvdwc(l,j)=gvdwc(l,j)+gg(l)
840 C now derivative over rik
841 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
842 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
847 gvdwx(m,i)=gvdwx(m,i)-gg(m)
848 gvdwx(m,k)=gvdwx(m,k)+gg(m)
851 gvdwc(l,i)=gvdwc(l,i)-gg(l)
852 gvdwc(l,k)=gvdwc(l,k)+gg(l)
854 C now derivative over rjk
855 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
856 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
861 gvdwx(m,j)=gvdwx(m,j)-gg(m)
862 gvdwx(m,k)=gvdwx(m,k)+gg(m)
865 gvdwc(l,j)=gvdwc(l,j)-gg(l)
866 gvdwc(l,k)=gvdwc(l,k)+gg(l)