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 'COMMON.SBRIDGE'
89 include 'COMMON.CHAIN'
90 include 'COMMON.DERIV'
91 include 'COMMON.LOCAL'
92 include 'COMMON.INTERACT'
94 include 'COMMON.IOUNITS'
96 include 'COMMON.NAMES'
104 double precision h_base
115 c integer itypi,itypj,k,l
116 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
117 double precision sig0ij,ljd,sig,fac,e1,e2
118 double precision dcosom1(3),dcosom2(3),ed
119 double precision pom1,pom2
120 double precision ljA,ljB,ljXs
121 double precision d_ljB(1:3)
122 double precision ssA,ssB,ssC,ssXs
123 double precision ssxm,ljxm,ssm,ljm
124 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
125 double precision f1,f2,h1,h2,hd1,hd2
126 double precision omega,delta_inv,deltasq_inv,fac1,fac2
128 double precision xm,d_xm(1:3)
129 double precision sslipi,sslipj,ssgradlipi,ssgradlipj
130 integer ici,icj,itypi,itypj
131 double precision boxshift,sscale,sscagrad
132 double precision aa,bb
133 c-------END FIRST METHOD
134 c-------SECOND METHOD
135 c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
136 c-------END SECOND METHOD
139 logical checkstop,transgrad
140 common /sschecks/ checkstop,transgrad
142 integer icheck,nicheck,jcheck,njcheck
143 double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
144 c-------END TESTING CODE
151 if (ici.eq.0 .or. icj.eq.0) then
152 write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
153 & "Attempt to create",
154 & " a disulfide link between non-cysteine residues ",restyp(i),i,
159 dxi=dc_norm(1,nres+i)
160 dyi=dc_norm(2,nres+i)
161 dzi=dc_norm(3,nres+i)
162 dsci_inv=vbld_inv(i+nres)
166 call to_box(xi,yi,zi)
167 C define scaling factor for lipids
169 C if (positi.le.0) positi=positi+boxzsize
171 C first for peptide groups
172 c for each residue check if it is in lipid or lipid water border area
173 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
178 call to_box(xj,yj,zj)
179 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
180 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
181 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
182 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
183 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
184 xj=boxshift(xj-xi,boxxsize)
185 yj=boxshift(yj-yi,boxysize)
186 zj=boxshift(zj-zi,boxzsize)
187 dxj=dc_norm(1,nres+j)
188 dyj=dc_norm(2,nres+j)
189 dzj=dc_norm(3,nres+j)
190 dscj_inv=vbld_inv(j+nres)
192 chi1=chi(itypi,itypj)
193 chi2=chi(itypj,itypi)
200 alf12=0.5D0*(alf1+alf2)
202 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
203 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
204 sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
205 sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
206 c The following are set in sc_angular
210 c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
211 c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
212 c om12=dxi*dxj+dyi*dyj+dzi*dzj
214 rij=1.0D0/rij ! Reset this so it makes sense
216 sig0ij=sigma(itypi,itypj)
217 sig=sig0ij*dsqrt(1.0D0/sigsq)
220 ljA=eps1*eps2rt**2*eps3rt**2
223 ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
228 deltat12=om2-om1+2.0d0
233 & +akth*(deltat1*deltat1+deltat2*deltat2)
234 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
235 ssxm=ssXs-0.5D0*ssB/ssA
238 c$$$c Some extra output
239 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
240 c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
241 c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
242 c$$$ if (ssx0.gt.0.0d0) then
243 c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
247 c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
248 c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
249 c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
251 c-------END TESTING CODE
254 c Stop and plot energy and derivative as a function of distance
256 ssm=ssC-0.25D0*ssB*ssB/ssA
257 ljm=-0.25D0*ljB*bb/aa
259 & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
267 if (.not.checkstop) then
274 if (checkstop) rij=(ssxm-1.0d0)+
275 & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
276 c-------END TESTING CODE
278 if (rij.gt.ljxm) then
281 fac=(1.0D0/ljd)**expon
284 eij=eps1*eps2rt*eps3rt*(e1+e2)
287 eij=eij*eps2rt*eps3rt*sss
290 e1=e1*eps1*eps2rt**2*eps3rt**2
291 ed=-expon*(e1+eij)/ljd
293 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
294 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
295 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
296 eom12=eij*eps1_om12+eps2der*eps2rt_om12
297 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
298 else if (rij.lt.ssxm) then
301 eij=ssA*ssd*ssd+ssB*ssd+ssC
303 ed=2*akcm*ssd+akct*deltat12
304 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
306 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
307 eom1=-2*akth*deltat1-pom1-om2*pom2
308 eom2= 2*akth*deltat2+pom1-om1*pom2
311 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
313 d_ssxm(1)=0.5D0*akct/ssA
317 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
318 d_ljxm(2)=d_ljxm(1)*sigsq_om2
319 d_ljxm(3)=d_ljxm(1)*sigsq_om12
320 d_ljxm(1)=d_ljxm(1)*sigsq_om1
322 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
325 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
329 ssm=ssC-0.25D0*ssB*ssB/ssA
330 d_ssm(1)=0.5D0*akct*ssB/ssA
331 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
332 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
334 f1=(rij-xm)/(ssxm-xm)
335 f2=(rij-ssxm)/(xm-ssxm)
339 delta_inv=1.0d0/(xm-ssxm)
340 deltasq_inv=delta_inv*delta_inv
342 fac1=deltasq_inv*fac*(xm-rij)
343 fac2=deltasq_inv*fac*(rij-ssxm)
344 ed=delta_inv*(Ht*hd2-ssm*hd1)
346 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
347 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
348 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
349 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
352 ljm=-0.25D0*ljB*bb/aa
353 d_ljm(1)=-0.5D0*bb/aa*ljB
354 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
355 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
357 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
358 f1=(rij-ljxm)/(xm-ljxm)
359 f2=(rij-xm)/(ljxm-xm)
363 delta_inv=1.0d0/(ljxm-xm)
364 deltasq_inv=delta_inv*delta_inv
366 fac1=deltasq_inv*fac*(ljxm-rij)
367 fac2=deltasq_inv*fac*(rij-xm)
368 ed=delta_inv*(ljm*hd2-Ht*hd1)
370 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
371 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
372 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
373 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
375 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
377 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
383 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
384 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
385 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
387 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
388 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
389 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
390 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
393 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
395 c$$$ d_ljm(k)=ljm*d_ljB(k)
399 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
400 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
401 c$$$ d_ss(2)=akct*ssd
402 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
403 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
406 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
407 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
408 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
410 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
411 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
413 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
415 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
416 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
417 c$$$ h1=h_base(f1,hd1)
418 c$$$ h2=h_base(f2,hd2)
419 c$$$ eij=ss*h1+ljf*h2
420 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
421 c$$$ deltasq_inv=delta_inv*delta_inv
422 c$$$ fac=ljf*hd2-ss*hd1
423 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
424 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
425 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
426 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
427 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
428 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
429 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
431 c$$$ havebond=.false.
432 c$$$ if (ed.gt.0.0d0) havebond=.true.
433 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
440 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
441 c write(iout,'(a15,f12.2,f8.1,2i5)')
442 c & "SSBOND_E_FORM",totT,t_bath,i,j
446 dyn_ssbond_ij(ici,icj)=eij
447 else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300)
449 dyn_ssbond_ij(ici,icj)=1.0d300
452 c write(iout,'(a15,f12.2,f8.1,2i5)')
453 c & "SSBOND_E_BREAK",totT,t_bath,i,j
460 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
461 & "CHECKSTOP",rij,eij,ed
466 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
473 c-------END TESTING CODE
474 gg_lipi(3)=ssgradlipi*eij
475 gg_lipj(3)=ssgradlipj*eij
478 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
479 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
482 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
485 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
486 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
487 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
488 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
489 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
490 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
494 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
499 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
500 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
505 C-----------------------------------------------------------------------------
507 double precision function h_base(x,deriv)
508 c A smooth function going 0->1 in range [0,1]
509 c It should NOT be called outside range [0,1], it will not work there.
516 double precision deriv
522 c Two parabolas put together. First derivative zero at extrema
523 c$$$ if (x.lt.0.5D0) then
524 c$$$ h_base=2.0D0*x*x
528 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
529 c$$$ deriv=4.0D0*deriv
532 c Third degree polynomial. First derivative zero at extrema
533 h_base=x*x*(3.0d0-2.0d0*x)
534 deriv=6.0d0*x*(1.0d0-x)
536 c Fifth degree polynomial. First and second derivatives zero at extrema
538 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
540 c$$$ deriv=deriv*deriv
541 c$$$ deriv=30.0d0*xsq*deriv
545 c----------------------------------------------------------------------------
546 subroutine dyn_set_nss
547 c Adjust nss and other relevant variables based on dyn_ssbond_ij
555 include 'COMMON.SBRIDGE'
556 include 'COMMON.CHAIN'
557 include 'COMMON.IOUNITS'
558 C include 'COMMON.SETUP'
561 C include 'COMMON.MD'
566 double precision emin
568 integer diff,allflag(maxdim_cont),allnss,
569 & allihpb(maxdim_cont),alljhpb(maxdim_cont),
570 & newnss,newihpb(maxdim_cont),newjhpb(maxdim_cont)
572 integer i_newnss(1024),displ(0:1024)
573 integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss
578 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
587 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
591 if (allflag(i).eq.0 .and.
592 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
593 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
597 if (emin.lt.1.0d300) then
600 if (allflag(i).eq.0 .and.
601 & (allihpb(i).eq.allihpb(imin) .or.
602 & alljhpb(i).eq.allihpb(imin) .or.
603 & allihpb(i).eq.alljhpb(imin) .or.
604 & alljhpb(i).eq.alljhpb(imin))) then
611 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
615 if (allflag(i).eq.1) then
617 newihpb(newnss)=allihpb(i)
618 newjhpb(newnss)=alljhpb(i)
623 if (nfgtasks.gt.1)then
625 call MPI_Reduce(newnss,g_newnss,1,
626 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
627 call MPI_Gather(newnss,1,MPI_INTEGER,
628 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
631 displ(i)=i_newnss(i-1)+displ(i-1)
633 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
634 & g_newihpb,i_newnss,displ,MPI_INTEGER,
636 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
637 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
639 if(fg_rank.eq.0) then
640 c print *,'g_newnss',g_newnss
641 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
642 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
645 newihpb(i)=g_newihpb(i)
646 newjhpb(i)=g_newjhpb(i)
654 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
659 if (idssb(i).eq.newihpb(j) .and.
660 & jdssb(i).eq.newjhpb(j)) found=.true.
664 c if (.not.found.and.fg_rank.eq.0)
665 c & write(iout,'(a15,f12.2,f8.1,2i5)')
666 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
674 if (newihpb(i).eq.idssb(j) .and.
675 & newjhpb(i).eq.jdssb(j)) found=.true.
679 c if (.not.found.and.fg_rank.eq.0)
680 c & write(iout,'(a15,f12.2,f8.1,2i5)')
681 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
695 c----------------------------------------------------------------------------
699 subroutine read_ssHist
704 include "DIMENSIONS.FREE"
705 include 'COMMON.FREE'
709 character*80 controlcard
712 call card_concat(controlcard,.true.)
714 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
721 c$$$C----------------------------------------------------------------------------
722 subroutine triple_ssbond_ene(resi,resj,resk,eij)
724 include 'COMMON.SBRIDGE'
725 include 'COMMON.CHAIN'
726 include 'COMMON.DERIV'
727 include 'COMMON.LOCAL'
728 include 'COMMON.INTERACT'
730 include 'COMMON.IOUNITS'
731 include 'COMMON.CALC'
734 C include 'COMMON.MD'
739 double precision h_base
743 integer resi,resj,resk
746 double precision eij,eij1,eij2,eij3
750 c integer itypi,itypj,k,l
751 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
752 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
753 double precision xik,yik,zik,xjk,yjk,zjk
754 double precision sig0ij,ljd,sig,fac,e1,e2
755 double precision dcosom1(3),dcosom2(3),ed
756 double precision pom1,pom2
757 double precision ljA,ljB,ljXs
758 double precision d_ljB(1:3)
759 double precision ssA,ssB,ssC,ssXs
760 double precision ssxm,ljxm,ssm,ljm
761 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
766 C write(iout,*) resi,resj,resk
768 dxi=dc_norm(1,nres+i)
769 dyi=dc_norm(2,nres+i)
770 dzi=dc_norm(3,nres+i)
771 dsci_inv=vbld_inv(i+nres)
781 dxj=dc_norm(1,nres+j)
782 dyj=dc_norm(2,nres+j)
783 dzj=dc_norm(3,nres+j)
784 dscj_inv=vbld_inv(j+nres)
790 dxk=dc_norm(1,nres+k)
791 dyk=dc_norm(2,nres+k)
792 dzk=dc_norm(3,nres+k)
793 dscj_inv=vbld_inv(k+nres)
803 rrij=(xij*xij+yij*yij+zij*zij)
804 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
805 rrik=(xik*xik+yik*yik+zik*zik)
807 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
809 C there are three combination of distances for each trisulfide bonds
810 C The first case the ith atom is the center
811 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
812 C distance y is second distance the a,b,c,d are parameters derived for
813 C this problem d parameter was set as a penalty currenlty set to 1.
814 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
815 C second case jth atom is center
816 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
817 C the third case kth atom is the center
818 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
823 C write(iout,*)i,j,k,eij
824 C The energy penalty calculated now time for the gradient part
825 C derivative over rij
826 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
827 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
832 gvdwx(m,i)=gvdwx(m,i)-gg(m)
833 gvdwx(m,j)=gvdwx(m,j)+gg(m)
836 gvdwc(l,i)=gvdwc(l,i)-gg(l)
837 gvdwc(l,j)=gvdwc(l,j)+gg(l)
839 C now derivative over rik
840 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
841 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
846 gvdwx(m,i)=gvdwx(m,i)-gg(m)
847 gvdwx(m,k)=gvdwx(m,k)+gg(m)
850 gvdwc(l,i)=gvdwc(l,i)-gg(l)
851 gvdwc(l,k)=gvdwc(l,k)+gg(l)
853 C now derivative over rjk
854 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
855 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
860 gvdwx(m,j)=gvdwx(m,j)-gg(m)
861 gvdwx(m,k)=gvdwx(m,k)+gg(m)
864 gvdwc(l,j)=gvdwc(l,j)-gg(l)
865 gvdwc(l,k)=gvdwc(l,k)+gg(l)