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 (iout,'(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 c write (iout,'(2(a,i5),4(a,f7.2))') "resi",resi," resj",resj,
280 c & " ljxm",ljxm," ljxs",ljxs," ssxm",ssxm," rij",rij
281 if (rij.gt.ljxm) then
284 fac=(1.0D0/ljd)**expon
287 eij=eps1*eps2rt*eps3rt*(e1+e2)
290 eij=eij*eps2rt*eps3rt*sss
293 e1=e1*eps1*eps2rt**2*eps3rt**2
294 ed=-expon*(e1+eij)/ljd
296 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
297 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
298 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
299 eom12=eij*eps1_om12+eps2der*eps2rt_om12
300 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
301 else if (rij.lt.ssxm) then
304 c write (iout,*) "ssMD: nssbond",nssbond
306 eij=ssA*ssd*ssd+ssB*ssd+ssC
308 ed=2*akcm*ssd+akct*deltat12
309 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
311 pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
312 eom1=-2*akth*deltat1-pom1-om2*pom2
313 eom2= 2*akth*deltat2+pom1-om1*pom2
317 omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
319 d_ssxm(1)=0.5D0*akct/ssA
323 d_ljxm(1)=sig0ij/sqrt(sigsq**3)
324 d_ljxm(2)=d_ljxm(1)*sigsq_om2
325 d_ljxm(3)=d_ljxm(1)*sigsq_om12
326 d_ljxm(1)=d_ljxm(1)*sigsq_om1
328 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
331 d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
335 ssm=ssC-0.25D0*ssB*ssB/ssA
336 d_ssm(1)=0.5D0*akct*ssB/ssA
337 d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
338 d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
340 f1=(rij-xm)/(ssxm-xm)
341 f2=(rij-ssxm)/(xm-ssxm)
345 delta_inv=1.0d0/(xm-ssxm)
346 deltasq_inv=delta_inv*delta_inv
348 fac1=deltasq_inv*fac*(xm-rij)
349 fac2=deltasq_inv*fac*(rij-ssxm)
350 ed=delta_inv*(Ht*hd2-ssm*hd1)
352 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
353 eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
354 eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
355 eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
358 ljm=-0.25D0*ljB*bb/aa
359 d_ljm(1)=-0.5D0*bb/aa*ljB
360 d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
361 d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
363 d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
364 f1=(rij-ljxm)/(xm-ljxm)
365 f2=(rij-xm)/(ljxm-xm)
369 delta_inv=1.0d0/(ljxm-xm)
370 deltasq_inv=delta_inv*delta_inv
372 fac1=deltasq_inv*fac*(ljxm-rij)
373 fac2=deltasq_inv*fac*(rij-xm)
374 ed=delta_inv*(ljm*hd2-Ht*hd1)
376 ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
377 eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
378 eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
379 eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
381 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
383 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
389 c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
390 c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
391 c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
393 c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
394 c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
395 c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
396 c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
399 c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
401 c$$$ d_ljm(k)=ljm*d_ljB(k)
405 c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
406 c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
407 c$$$ d_ss(2)=akct*ssd
408 c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
409 c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
412 c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
413 c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
414 c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
416 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
417 c$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
419 c$$$ ljf=ljm+ljf*ljB*fac1*fac1
421 c$$$ f1=(rij-ljxm)/(ssxm-ljxm)
422 c$$$ f2=(rij-ssxm)/(ljxm-ssxm)
423 c$$$ h1=h_base(f1,hd1)
424 c$$$ h2=h_base(f2,hd2)
425 c$$$ eij=ss*h1+ljf*h2
426 c$$$ delta_inv=1.0d0/(ljxm-ssxm)
427 c$$$ deltasq_inv=delta_inv*delta_inv
428 c$$$ fac=ljf*hd2-ss*hd1
429 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
430 c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
431 c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
432 c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
433 c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
434 c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
435 c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
437 c$$$ havebond=.false.
438 c$$$ if (ed.gt.0.0d0) havebond=.true.
439 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
446 c if (dyn_ssbond_ij(i,j).eq.1.0d300) then
447 c write(iout,'(a15,f12.2,f8.1,2i5)')
448 c & "SSBOND_E_FORM",totT,t_bath,i,j
452 dyn_ssbond_ij(ici,icj)=eij
453 else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300)
455 dyn_ssbond_ij(ici,icj)=1.0d300
458 c write(iout,'(a15,f12.2,f8.1,2i5)')
459 c & "SSBOND_E_BREAK",totT,t_bath,i,j
466 if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
467 & "CHECKSTOP",rij,eij,ed
472 write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
479 c-------END TESTING CODE
480 gg_lipi(3)=ssgradlipi*eij
481 gg_lipj(3)=ssgradlipj*eij
484 dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
485 dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
488 gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
491 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
492 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
493 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
494 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
495 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
496 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
500 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
505 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
506 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
511 C-----------------------------------------------------------------------------
513 double precision function h_base(x,deriv)
514 c A smooth function going 0->1 in range [0,1]
515 c It should NOT be called outside range [0,1], it will not work there.
522 double precision deriv
528 c Two parabolas put together. First derivative zero at extrema
529 c$$$ if (x.lt.0.5D0) then
530 c$$$ h_base=2.0D0*x*x
534 c$$$ h_base=1.0D0-2.0D0*deriv*deriv
535 c$$$ deriv=4.0D0*deriv
538 c Third degree polynomial. First derivative zero at extrema
539 h_base=x*x*(3.0d0-2.0d0*x)
540 deriv=6.0d0*x*(1.0d0-x)
542 c Fifth degree polynomial. First and second derivatives zero at extrema
544 c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
546 c$$$ deriv=deriv*deriv
547 c$$$ deriv=30.0d0*xsq*deriv
551 c----------------------------------------------------------------------------
552 subroutine dyn_set_nss
553 c Adjust nss and other relevant variables based on dyn_ssbond_ij
561 include 'COMMON.SBRIDGE'
562 include 'COMMON.CHAIN'
563 include 'COMMON.IOUNITS'
566 C include 'COMMON.MD'
571 double precision emin
573 integer diff,allflag(maxdim_cont),allnss,
574 & allihpb(maxdim_cont),alljhpb(maxdim_cont),
575 & newnss,newihpb(maxdim_cont),newjhpb(maxdim_cont)
577 integer i_newnss(1024),displ(0:1024)
578 integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss
584 if (dyn_ssbond_ij(i,j).lt.1.0d300) then
593 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
597 if (allflag(i).eq.0 .and.
598 & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
599 emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
603 if (emin.lt.1.0d300) then
606 if (allflag(i).eq.0 .and.
607 & (allihpb(i).eq.allihpb(imin) .or.
608 & alljhpb(i).eq.allihpb(imin) .or.
609 & allihpb(i).eq.alljhpb(imin) .or.
610 & alljhpb(i).eq.alljhpb(imin))) then
617 cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
621 if (allflag(i).eq.1) then
623 newihpb(newnss)=allihpb(i)
624 newjhpb(newnss)=alljhpb(i)
629 if (nfgtasks.gt.1)then
631 call MPI_Reduce(newnss,g_newnss,1,
632 & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
633 call MPI_Gather(newnss,1,MPI_INTEGER,
634 & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
637 displ(i)=i_newnss(i-1)+displ(i-1)
639 call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
640 & g_newihpb,i_newnss,displ,MPI_INTEGER,
642 call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
643 & g_newjhpb,i_newnss,displ,MPI_INTEGER,
645 if(fg_rank.eq.0) then
646 c print *,'g_newnss',g_newnss
647 c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
648 c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
651 newihpb(i)=g_newihpb(i)
652 newjhpb(i)=g_newjhpb(i)
660 cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
665 if (idssb(i).eq.newihpb(j) .and.
666 & jdssb(i).eq.newjhpb(j)) found=.true.
670 c if (.not.found.and.fg_rank.eq.0)
671 c & write(iout,'(a15,f12.2,f8.1,2i5)')
672 c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
680 if (newihpb(i).eq.idssb(j) .and.
681 & newjhpb(i).eq.jdssb(j)) found=.true.
685 c if (.not.found.and.fg_rank.eq.0)
686 c & write(iout,'(a15,f12.2,f8.1,2i5)')
687 c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
701 c----------------------------------------------------------------------------
705 subroutine read_ssHist
710 include "DIMENSIONS.FREE"
711 include 'COMMON.FREE'
715 character*80 controlcard
718 call card_concat(controlcard,.true.)
720 & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
727 c$$$C----------------------------------------------------------------------------
728 subroutine triple_ssbond_ene(resi,resj,resk,eij)
730 include 'COMMON.SBRIDGE'
731 include 'COMMON.CHAIN'
732 include 'COMMON.DERIV'
733 include 'COMMON.LOCAL'
734 include 'COMMON.INTERACT'
736 include 'COMMON.IOUNITS'
737 include 'COMMON.CALC'
740 C include 'COMMON.MD'
745 double precision h_base
749 integer resi,resj,resk
752 double precision eij,eij1,eij2,eij3
756 c integer itypi,itypj,k,l
757 double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
758 double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
759 double precision xik,yik,zik,xjk,yjk,zjk
760 double precision sig0ij,ljd,sig,fac,e1,e2
761 double precision dcosom1(3),dcosom2(3),ed
762 double precision pom1,pom2
763 double precision ljA,ljB,ljXs
764 double precision d_ljB(1:3)
765 double precision ssA,ssB,ssC,ssXs
766 double precision ssxm,ljxm,ssm,ljm
767 double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
772 C write(iout,*) resi,resj,resk
774 dxi=dc_norm(1,nres+i)
775 dyi=dc_norm(2,nres+i)
776 dzi=dc_norm(3,nres+i)
777 dsci_inv=vbld_inv(i+nres)
787 dxj=dc_norm(1,nres+j)
788 dyj=dc_norm(2,nres+j)
789 dzj=dc_norm(3,nres+j)
790 dscj_inv=vbld_inv(j+nres)
796 dxk=dc_norm(1,nres+k)
797 dyk=dc_norm(2,nres+k)
798 dzk=dc_norm(3,nres+k)
799 dscj_inv=vbld_inv(k+nres)
809 rrij=(xij*xij+yij*yij+zij*zij)
810 rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
811 rrik=(xik*xik+yik*yik+zik*zik)
813 rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
815 C there are three combination of distances for each trisulfide bonds
816 C The first case the ith atom is the center
817 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
818 C distance y is second distance the a,b,c,d are parameters derived for
819 C this problem d parameter was set as a penalty currenlty set to 1.
820 eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
821 C second case jth atom is center
822 eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
823 C the third case kth atom is the center
824 eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
829 C write(iout,*)i,j,k,eij
830 C The energy penalty calculated now time for the gradient part
831 C derivative over rij
832 fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
833 &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))
838 gvdwx(m,i)=gvdwx(m,i)-gg(m)
839 gvdwx(m,j)=gvdwx(m,j)+gg(m)
842 gvdwc(l,i)=gvdwc(l,i)-gg(l)
843 gvdwc(l,j)=gvdwc(l,j)+gg(l)
845 C now derivative over rik
846 fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
847 &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
852 gvdwx(m,i)=gvdwx(m,i)-gg(m)
853 gvdwx(m,k)=gvdwx(m,k)+gg(m)
856 gvdwc(l,i)=gvdwc(l,i)-gg(l)
857 gvdwc(l,k)=gvdwc(l,k)+gg(l)
859 C now derivative over rjk
860 fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
861 &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
866 gvdwx(m,j)=gvdwx(m,j)-gg(m)
867 gvdwx(m,k)=gvdwx(m,k)+gg(m)
870 gvdwc(l,j)=gvdwc(l,j)-gg(l)
871 gvdwc(l,k)=gvdwc(l,k)+gg(l)