ctest checkgrad ref ene and limits
[unres.git] / source / cluster / wham / src-M / ssMD.F
1 c----------------------------------------------------------------------------
2       subroutine check_energies
3 c      implicit none
4
5 c     Includes
6       include 'DIMENSIONS'
7       include 'COMMON.CHAIN'
8       include 'COMMON.VAR'
9       include 'COMMON.IOUNITS'
10       include 'COMMON.SBRIDGE'
11       include 'COMMON.LOCAL'
12       include 'COMMON.GEO'
13
14 c     External functions
15       double precision ran_number
16       external ran_number
17
18 c     Local variables
19       integer i,j,k,l,lmax,p,pmax
20       double precision rmin,rmax
21       double precision eij
22
23       double precision d
24       double precision wi,rij,tj,pj
25
26
27 c      return
28
29       i=5
30       j=14
31
32       d=dsc(1)
33       rmin=2.0D0
34       rmax=12.0D0
35
36       lmax=10000
37       pmax=1
38
39       do k=1,3
40         c(k,i)=0.0D0
41         c(k,j)=0.0D0
42         c(k,nres+i)=0.0D0
43         c(k,nres+j)=0.0D0
44       enddo
45
46       do l=1,lmax
47
48 ct        wi=ran_number(0.0D0,pi)
49 c        wi=ran_number(0.0D0,pi/6.0D0)
50 c        wi=0.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)
54 c        pj=0.0D0
55
56         do p=1,pmax
57 ct           rij=ran_number(rmin,rmax)
58
59            c(1,j)=d*sin(pj)*cos(tj)
60            c(2,j)=d*sin(pj)*sin(tj)
61            c(3,j)=d*cos(pj)
62
63            c(3,nres+i)=-rij
64
65            c(1,i)=d*sin(wi)
66            c(3,i)=-rij-d*cos(wi)
67
68            do k=1,3
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
73            enddo
74
75            call dyn_ssbond_ene(i,j,eij)
76         enddo
77       enddo
78
79       call exit(1)
80
81       return
82       end
83
84 C-----------------------------------------------------------------------------
85
86       subroutine dyn_ssbond_ene(resi,resj,eij)
87 c      implicit none
88
89 c     Includes
90       include 'DIMENSIONS'
91       include 'COMMON.SBRIDGE'
92       include 'COMMON.CHAIN'
93       include 'COMMON.DERIV'
94       include 'COMMON.LOCAL'
95       include 'COMMON.INTERACT'
96       include 'COMMON.VAR'
97       include 'COMMON.IOUNITS'
98       include 'COMMON.CALC'
99 #ifndef CLUST
100 #ifndef WHAM
101 C      include 'COMMON.MD'
102 #endif
103 #endif
104
105 c     External functions
106       double precision h_base
107       external h_base
108
109 c     Input arguments
110       integer resi,resj
111
112 c     Output arguments
113       double precision eij
114
115 c     Local variables
116       logical havebond
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
129 c-------FIRST METHOD
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
135
136 c-------TESTING CODE
137       logical checkstop,transgrad
138       common /sschecks/ checkstop,transgrad
139
140       integer icheck,nicheck,jcheck,njcheck
141       double precision echeck(-1:1),deps,ssx0,ljx0
142 c-------END TESTING CODE
143
144
145       i=resi
146       j=resj
147
148       itypi=itype(i)
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)
153         xi=c(1,nres+i)
154         yi=c(2,nres+i)
155         zi=c(3,nres+i)
156           xi=dmod(xi,boxxsize)
157           if (xi.lt.0) xi=xi+boxxsize
158           yi=dmod(yi,boxysize)
159           if (yi.lt.0) yi=yi+boxysize
160           zi=dmod(zi,boxzsize)
161           if (zi.lt.0) zi=zi+boxzsize
162 C define scaling factor for lipids
163
164 C        if (positi.le.0) positi=positi+boxzsize
165 C        print *,i
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
173          fracinbuf=1.0d0-
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
182         else
183          sslipi=1.0d0
184          ssgradlipi=0.0
185         endif
186        else
187          sslipi=0.0d0
188          ssgradlipi=0.0
189        endif
190
191       itypj=itype(j)
192             xj=c(1,nres+j)
193             yj=c(2,nres+j)
194             zj=c(3,nres+j)
195           xj=dmod(xj,boxxsize)
196           if (xj.lt.0) xj=xj+boxxsize
197           yj=dmod(yj,boxysize)
198           if (yj.lt.0) yj=yj+boxysize
199           zj=dmod(zj,boxzsize)
200           if (zj.lt.0) zj=zj+boxzsize
201        if ((zj.gt.bordlipbot)
202      &.and.(zj.lt.bordliptop)) then
203 C the energy transfer exist
204         if (zj.lt.buflipbot) then
205 C what fraction I am in
206          fracinbuf=1.0d0-
207      &        ((positi-bordlipbot)/lipbufthick)
208 C lipbufthick is thickenes of lipid buffore
209          sslipj=sscalelip(fracinbuf)
210          ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
211         elseif (zi.gt.bufliptop) then
212          fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
213          sslipj=sscalelip(fracinbuf)
214          ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
215         else
216          sslipj=1.0d0
217          ssgradlipj=0.0
218         endif
219        else
220          sslipj=0.0d0
221          ssgradlipj=0.0
222        endif
223       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
224      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
225       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
226      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
227
228       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
229       xj_safe=xj
230       yj_safe=yj
231       zj_safe=zj
232       subchap=0
233       xj_safe=xj
234       yj_safe=yj
235       zj_safe=zj
236       subchap=0
237       do xshift=-1,1
238       do yshift=-1,1
239       do zshift=-1,1
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           if(dist_temp.lt.dist_init) then
245             dist_init=dist_temp
246             xj_temp=xj
247             yj_temp=yj
248             zj_temp=zj
249             subchap=1
250           endif
251        enddo
252        enddo
253        enddo
254        if (subchap.eq.1) then
255           xj=xj_temp-xi
256           yj=yj_temp-yi
257           zj=zj_temp-zi
258        else
259           xj=xj_safe-xi
260           yj=yj_safe-yi
261           zj=zj_safe-zi
262        endif
263
264       dxj=dc_norm(1,nres+j)
265       dyj=dc_norm(2,nres+j)
266       dzj=dc_norm(3,nres+j)
267       dscj_inv=vbld_inv(j+nres)
268
269       chi1=chi(itypi,itypj)
270       chi2=chi(itypj,itypi)
271       chi12=chi1*chi2
272       chip1=chip(itypi)
273       chip2=chip(itypj)
274       chip12=chip1*chip2
275       alf1=alp(itypi)
276       alf2=alp(itypj)
277       alf12=0.5D0*(alf1+alf2)
278
279       rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
280       rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
281 c     The following are set in sc_angular
282 c      erij(1)=xj*rij
283 c      erij(2)=yj*rij
284 c      erij(3)=zj*rij
285 c      om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
286 c      om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
287 c      om12=dxi*dxj+dyi*dyj+dzi*dzj
288       call sc_angular
289       rij=1.0D0/rij  ! Reset this so it makes sense
290
291       sig0ij=sigma(itypi,itypj)
292       sig=sig0ij*dsqrt(1.0D0/sigsq)
293
294       ljXs=sig-sig0ij
295       ljA=eps1*eps2rt**2*eps3rt**2
296       ljB=ljA*bb
297       ljA=ljA*aa
298       ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
299
300       ssXs=d0cm
301       deltat1=1.0d0-om1
302       deltat2=1.0d0+om2
303       deltat12=om2-om1+2.0d0
304       cosphi=om12-om1*om2
305       ssA=akcm
306       ssB=akct*deltat12
307       ssC=ss_depth
308      &     +akth*(deltat1*deltat1+deltat2*deltat2)
309      &     +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
310       ssxm=ssXs-0.5D0*ssB/ssA
311
312 c-------TESTING CODE
313 c$$$c     Some extra output
314 c$$$      ssm=ssC-0.25D0*ssB*ssB/ssA
315 c$$$      ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
316 c$$$      ssx0=ssB*ssB-4.0d0*ssA*ssC
317 c$$$      if (ssx0.gt.0.0d0) then
318 c$$$        ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
319 c$$$      else
320 c$$$        ssx0=ssxm
321 c$$$      endif
322 c$$$      ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
323 c$$$      write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
324 c$$$     &     ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
325 c$$$      return
326 c-------END TESTING CODE
327
328 c-------TESTING CODE
329 c     Stop and plot energy and derivative as a function of distance
330       if (checkstop) then
331         ssm=ssC-0.25D0*ssB*ssB/ssA
332         ljm=-0.25D0*ljB*bb/aa
333         if (ssm.lt.ljm .and.
334      &       dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
335           nicheck=1000
336           njcheck=1
337           deps=0.5d-7
338         else
339           checkstop=.false.
340         endif
341       endif
342       if (.not.checkstop) then
343         nicheck=0
344         njcheck=-1
345       endif
346
347       do icheck=0,nicheck
348       do jcheck=-1,njcheck
349       if (checkstop) rij=(ssxm-1.0d0)+
350      &       ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
351 c-------END TESTING CODE
352
353       if (rij.gt.ljxm) then
354         havebond=.false.
355         ljd=rij-ljXs
356         fac=(1.0D0/ljd)**expon
357         e1=fac*fac*aa
358         e2=fac*bb
359         eij=eps1*eps2rt*eps3rt*(e1+e2)
360 C        write(iout,*) eij,'TU?1'
361         eps2der=eij*eps3rt
362         eps3der=eij*eps2rt
363         eij=eij*eps2rt*eps3rt
364
365         sigder=-sig/sigsq
366         e1=e1*eps1*eps2rt**2*eps3rt**2
367         ed=-expon*(e1+eij)/ljd
368         sigder=ed*sigder
369         eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
370         eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
371         eom12=eij*eps1_om12+eps2der*eps2rt_om12
372      &       -2.0D0*alf12*eps3der+sigder*sigsq_om12
373       else if (rij.lt.ssxm) then
374         havebond=.true.
375         ssd=rij-ssXs
376         eij=ssA*ssd*ssd+ssB*ssd+ssC
377 C        write(iout,*) 'TU?2',ssc,ssd
378         ed=2*akcm*ssd+akct*deltat12
379         pom1=akct*ssd
380         pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
381         eom1=-2*akth*deltat1-pom1-om2*pom2
382         eom2= 2*akth*deltat2+pom1-om1*pom2
383         eom12=pom2
384       else
385         omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
386
387         d_ssxm(1)=0.5D0*akct/ssA
388         d_ssxm(2)=-d_ssxm(1)
389         d_ssxm(3)=0.0D0
390
391         d_ljxm(1)=sig0ij/sqrt(sigsq**3)
392         d_ljxm(2)=d_ljxm(1)*sigsq_om2
393         d_ljxm(3)=d_ljxm(1)*sigsq_om12
394         d_ljxm(1)=d_ljxm(1)*sigsq_om1
395
396 c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
397         xm=0.5d0*(ssxm+ljxm)
398         do k=1,3
399           d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
400         enddo
401         if (rij.lt.xm) then
402           havebond=.true.
403           ssm=ssC-0.25D0*ssB*ssB/ssA
404           d_ssm(1)=0.5D0*akct*ssB/ssA
405           d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
406           d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
407           d_ssm(3)=omega
408           f1=(rij-xm)/(ssxm-xm)
409           f2=(rij-ssxm)/(xm-ssxm)
410           h1=h_base(f1,hd1)
411           h2=h_base(f2,hd2)
412           eij=ssm*h1+Ht*h2
413 C         write(iout,*) eij,'TU?3'
414           delta_inv=1.0d0/(xm-ssxm)
415           deltasq_inv=delta_inv*delta_inv
416           fac=ssm*hd1-Ht*hd2
417           fac1=deltasq_inv*fac*(xm-rij)
418           fac2=deltasq_inv*fac*(rij-ssxm)
419           ed=delta_inv*(Ht*hd2-ssm*hd1)
420           eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
421           eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
422           eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
423         else
424           havebond=.false.
425           ljm=-0.25D0*ljB*bb/aa
426           d_ljm(1)=-0.5D0*bb/aa*ljB
427           d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
428           d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
429      +         alf12/eps3rt)
430           d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
431           f1=(rij-ljxm)/(xm-ljxm)
432           f2=(rij-xm)/(ljxm-xm)
433           h1=h_base(f1,hd1)
434           h2=h_base(f2,hd2)
435           eij=Ht*h1+ljm*h2
436 C        write(iout,*) 'TU?4',ssA
437           delta_inv=1.0d0/(ljxm-xm)
438           deltasq_inv=delta_inv*delta_inv
439           fac=Ht*hd1-ljm*hd2
440           fac1=deltasq_inv*fac*(ljxm-rij)
441           fac2=deltasq_inv*fac*(rij-xm)
442           ed=delta_inv*(ljm*hd2-Ht*hd1)
443           eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
444           eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
445           eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
446         endif
447 c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
448
449 c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
450 c$$$        ssd=rij-ssXs
451 c$$$        ljd=rij-ljXs
452 c$$$        fac1=rij-ljxm
453 c$$$        fac2=rij-ssxm
454 c$$$
455 c$$$        d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
456 c$$$        d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
457 c$$$        d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
458 c$$$
459 c$$$        ssm=ssC-0.25D0*ssB*ssB/ssA
460 c$$$        d_ssm(1)=0.5D0*akct*ssB/ssA
461 c$$$        d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
462 c$$$        d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
463 c$$$        d_ssm(3)=omega
464 c$$$
465 c$$$        ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
466 c$$$        do k=1,3
467 c$$$          d_ljm(k)=ljm*d_ljB(k)
468 c$$$        enddo
469 c$$$        ljm=ljm*ljB
470 c$$$
471 c$$$        ss=ssA*ssd*ssd+ssB*ssd+ssC
472 c$$$        d_ss(0)=2.0d0*ssA*ssd+ssB
473 c$$$        d_ss(2)=akct*ssd
474 c$$$        d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
475 c$$$        d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
476 c$$$        d_ss(3)=omega
477 c$$$
478 c$$$        ljf=bb(itypi,itypj)/aa(itypi,itypj)
479 c$$$        ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
480 c$$$        d_ljf(0)=ljf*2.0d0*ljB*fac1
481 c$$$        do k=1,3
482 c$$$          d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
483 c$$$     &         2.0d0*ljB*fac1*d_ljxm(k))
484 c$$$        enddo
485 c$$$        ljf=ljm+ljf*ljB*fac1*fac1
486 c$$$
487 c$$$        f1=(rij-ljxm)/(ssxm-ljxm)
488 c$$$        f2=(rij-ssxm)/(ljxm-ssxm)
489 c$$$        h1=h_base(f1,hd1)
490 c$$$        h2=h_base(f2,hd2)
491 c$$$        eij=ss*h1+ljf*h2
492 c$$$        delta_inv=1.0d0/(ljxm-ssxm)
493 c$$$        deltasq_inv=delta_inv*delta_inv
494 c$$$        fac=ljf*hd2-ss*hd1
495 c$$$        ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
496 c$$$        eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
497 c$$$     &       (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
498 c$$$        eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
499 c$$$     &       (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
500 c$$$        eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
501 c$$$     &       (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
502 c$$$
503 c$$$        havebond=.false.
504 c$$$        if (ed.gt.0.0d0) havebond=.true.
505 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
506
507       endif
508 C      write(iout,*) 'havebond',havebond
509       if (havebond) then
510 #ifndef CLUST
511 #ifndef WHAM
512 c        if (dyn_ssbond_ij(i,j).eq.1.0d300) then
513 c          write(iout,'(a15,f12.2,f8.1,2i5)')
514 c     &         "SSBOND_E_FORM",totT,t_bath,i,j
515 c        endif
516 #endif
517 #endif
518         dyn_ssbond_ij(i,j)=eij
519       else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
520         dyn_ssbond_ij(i,j)=1.0d300
521 #ifndef CLUST
522 #ifndef WHAM
523 c        write(iout,'(a15,f12.2,f8.1,2i5)')
524 c     &       "SSBOND_E_BREAK",totT,t_bath,i,j
525 #endif
526 #endif
527       endif
528
529 c-------TESTING CODE
530       if (checkstop) then
531         if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
532      &       "CHECKSTOP",rij,eij,ed
533         echeck(jcheck)=eij
534       endif
535       enddo
536       if (checkstop) then
537         write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
538       endif
539       enddo
540       if (checkstop) then
541         transgrad=.true.
542         checkstop=.false.
543       endif
544 c-------END TESTING CODE
545
546       do k=1,3
547         dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
548         dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
549       enddo
550       do k=1,3
551         gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
552       enddo
553       do k=1,3
554         gvdwx(k,i)=gvdwx(k,i)-gg(k)
555      &       +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
556      &       +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
557         gvdwx(k,j)=gvdwx(k,j)+gg(k)
558      &       +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
559      &       +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
560       enddo
561 cgrad      do k=i,j-1
562 cgrad        do l=1,3
563 cgrad          gvdwc(l,k)=gvdwc(l,k)+gg(l)
564 cgrad        enddo
565 cgrad      enddo
566
567       do l=1,3
568         gvdwc(l,i)=gvdwc(l,i)-gg(l)
569         gvdwc(l,j)=gvdwc(l,j)+gg(l)
570       enddo
571
572       return
573       end
574
575 C-----------------------------------------------------------------------------
576
577       double precision function h_base(x,deriv)
578 c     A smooth function going 0->1 in range [0,1]
579 c     It should NOT be called outside range [0,1], it will not work there.
580       implicit none
581
582 c     Input arguments
583       double precision x
584
585 c     Output arguments
586       double precision deriv
587
588 c     Local variables
589       double precision xsq
590
591
592 c     Two parabolas put together.  First derivative zero at extrema
593 c$$$      if (x.lt.0.5D0) then
594 c$$$        h_base=2.0D0*x*x
595 c$$$        deriv=4.0D0*x
596 c$$$      else
597 c$$$        deriv=1.0D0-x
598 c$$$        h_base=1.0D0-2.0D0*deriv*deriv
599 c$$$        deriv=4.0D0*deriv
600 c$$$      endif
601
602 c     Third degree polynomial.  First derivative zero at extrema
603       h_base=x*x*(3.0d0-2.0d0*x)
604       deriv=6.0d0*x*(1.0d0-x)
605
606 c     Fifth degree polynomial.  First and second derivatives zero at extrema
607 c$$$      xsq=x*x
608 c$$$      h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
609 c$$$      deriv=x-1.0d0
610 c$$$      deriv=deriv*deriv
611 c$$$      deriv=30.0d0*xsq*deriv
612
613       return
614       end
615
616 c----------------------------------------------------------------------------
617
618       subroutine dyn_set_nss
619 c     Adjust nss and other relevant variables based on dyn_ssbond_ij
620 c      implicit none
621
622 c     Includes
623       include 'DIMENSIONS'
624 #ifdef MPI
625       include "mpif.h"
626 #endif
627       include 'COMMON.SBRIDGE'
628       include 'COMMON.CHAIN'
629       include 'COMMON.IOUNITS'
630 C      include 'COMMON.SETUP'
631 #ifndef CLUST
632 #ifndef WHAM
633 C      include 'COMMON.MD'
634 #endif
635 #endif
636
637 c     Local variables
638       double precision emin
639       integer i,j,imin
640       integer diff,allflag(maxdim),allnss,
641      &     allihpb(maxdim),alljhpb(maxdim),
642      &     newnss,newihpb(maxdim),newjhpb(maxdim)
643       logical found
644       integer i_newnss(1024),displ(0:1024)
645       integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
646
647       allnss=0
648       do i=1,nres-1
649         do j=i+1,nres
650           if (dyn_ssbond_ij(i,j).lt.1.0d300) then
651             allnss=allnss+1
652             allflag(allnss)=0
653             allihpb(allnss)=i
654             alljhpb(allnss)=j
655           endif
656         enddo
657       enddo
658
659 cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
660
661  1    emin=1.0d300
662       do i=1,allnss
663         if (allflag(i).eq.0 .and.
664      &       dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
665           emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
666           imin=i
667         endif
668       enddo
669       if (emin.lt.1.0d300) then
670         allflag(imin)=1
671         do i=1,allnss
672           if (allflag(i).eq.0 .and.
673      &         (allihpb(i).eq.allihpb(imin) .or.
674      &         alljhpb(i).eq.allihpb(imin) .or.
675      &         allihpb(i).eq.alljhpb(imin) .or.
676      &         alljhpb(i).eq.alljhpb(imin))) then
677             allflag(i)=-1
678           endif
679         enddo
680         goto 1
681       endif
682
683 cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
684
685       newnss=0
686       do i=1,allnss
687         if (allflag(i).eq.1) then
688           newnss=newnss+1
689           newihpb(newnss)=allihpb(i)
690           newjhpb(newnss)=alljhpb(i)
691         endif
692       enddo
693
694 #ifdef MPI
695       if (nfgtasks.gt.1)then
696
697         call MPI_Reduce(newnss,g_newnss,1,
698      &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
699         call MPI_Gather(newnss,1,MPI_INTEGER,
700      &                  i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
701         displ(0)=0
702         do i=1,nfgtasks-1,1
703           displ(i)=i_newnss(i-1)+displ(i-1)
704         enddo
705         call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
706      &                   g_newihpb,i_newnss,displ,MPI_INTEGER,
707      &                   king,FG_COMM,IERR)     
708         call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
709      &                   g_newjhpb,i_newnss,displ,MPI_INTEGER,
710      &                   king,FG_COMM,IERR)     
711         if(fg_rank.eq.0) then
712 c         print *,'g_newnss',g_newnss
713 c         print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
714 c         print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
715          newnss=g_newnss  
716          do i=1,newnss
717           newihpb(i)=g_newihpb(i)
718           newjhpb(i)=g_newjhpb(i)
719          enddo
720         endif
721       endif
722 #endif
723
724       diff=newnss-nss
725
726 cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
727
728       do i=1,nss
729         found=.false.
730         do j=1,newnss
731           if (idssb(i).eq.newihpb(j) .and.
732      &         jdssb(i).eq.newjhpb(j)) found=.true.
733         enddo
734 #ifndef CLUST
735 #ifndef WHAM
736 c        if (.not.found.and.fg_rank.eq.0) 
737 c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
738 c     &       "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
739 #endif
740 #endif
741       enddo
742
743       do i=1,newnss
744         found=.false.
745         do j=1,nss
746           if (newihpb(i).eq.idssb(j) .and.
747      &         newjhpb(i).eq.jdssb(j)) found=.true.
748         enddo
749 #ifndef CLUST
750 #ifndef WHAM
751 c        if (.not.found.and.fg_rank.eq.0) 
752 c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
753 c     &       "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
754 #endif
755 #endif
756       enddo
757
758       nss=newnss
759       do i=1,nss
760         idssb(i)=newihpb(i)
761         jdssb(i)=newjhpb(i)
762       enddo
763
764       return
765       end
766
767
768 c$$$c-----------------------------------------------------------------------------
769 c$$$
770 c$$$      subroutine ss_relax(i_in,j_in)
771 c$$$      implicit none
772 c$$$
773 c$$$c     Includes
774 c$$$      include 'DIMENSIONS'
775 c$$$      include 'COMMON.VAR'
776 c$$$      include 'COMMON.CHAIN'
777 c$$$      include 'COMMON.IOUNITS'
778 c$$$      include 'COMMON.INTERACT'
779 c$$$
780 c$$$c     Input arguments
781 c$$$      integer i_in,j_in
782 c$$$
783 c$$$c     Local variables
784 c$$$      integer i,iretcode,nfun_sc
785 c$$$      logical scfail
786 c$$$      double precision var(maxvar),e_sc,etot
787 c$$$
788 c$$$
789 c$$$      mask_r=.true.
790 c$$$      do i=nnt,nct
791 c$$$        mask_side(i)=0
792 c$$$      enddo
793 c$$$      mask_side(i_in)=1
794 c$$$      mask_side(j_in)=1
795 c$$$
796 c$$$c     Minimize the two selected side-chains
797 c$$$      call overlap_sc(scfail)  ! Better not fail!
798 c$$$      call minimize_sc(e_sc,var,iretcode,nfun_sc)
799 c$$$
800 c$$$      mask_r=.false.
801 c$$$
802 c$$$      return
803 c$$$      end
804 c$$$
805 c$$$c-------------------------------------------------------------
806 c$$$
807 c$$$      subroutine minimize_sc(etot_sc,iretcode,nfun)
808 c$$$c     Minimize side-chains only, starting from geom but without modifying
809 c$$$c     bond lengths.
810 c$$$c     If mask_r is already set, only the selected side-chains are minimized,
811 c$$$c     otherwise all side-chains are minimized keeping the backbone frozen.
812 c$$$      implicit none
813 c$$$
814 c$$$c     Includes
815 c$$$      include 'DIMENSIONS'
816 c$$$      include 'COMMON.IOUNITS'
817 c$$$      include 'COMMON.VAR'
818 c$$$      include 'COMMON.CHAIN'
819 c$$$      include 'COMMON.GEO'
820 c$$$      include 'COMMON.MINIM'
821 c$$$      integer icall
822 c$$$      common /srutu/ icall
823 c$$$
824 c$$$c     Output arguments
825 c$$$      double precision etot_sc
826 c$$$      integer iretcode,nfun
827 c$$$
828 c$$$c     External functions/subroutines
829 c$$$      external func_sc,grad_sc,fdum
830 c$$$
831 c$$$c     Local variables
832 c$$$      integer liv,lv
833 c$$$      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
834 c$$$      integer iv(liv)
835 c$$$      double precision rdum(1)
836 c$$$      double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
837 c$$$      integer idum(1)
838 c$$$      integer i,nvar_restr
839 c$$$
840 c$$$
841 c$$$cmc      start_minim=.true.
842 c$$$      call deflt(2,iv,liv,lv,v)                                         
843 c$$$* 12 means fresh start, dont call deflt                                 
844 c$$$      iv(1)=12                                                          
845 c$$$* max num of fun calls                                                  
846 c$$$      if (maxfun.eq.0) maxfun=500
847 c$$$      iv(17)=maxfun
848 c$$$* max num of iterations                                                 
849 c$$$      if (maxmin.eq.0) maxmin=1000
850 c$$$      iv(18)=maxmin
851 c$$$* controls output                                                       
852 c$$$      iv(19)=1
853 c$$$* selects output unit                                                   
854 c$$$      iv(21)=0
855 c$$$c      iv(21)=iout               ! DEBUG
856 c$$$c      iv(21)=8                  ! DEBUG
857 c$$$* 1 means to print out result                                           
858 c$$$      iv(22)=0
859 c$$$c      iv(22)=1                  ! DEBUG
860 c$$$* 1 means to print out summary stats                                    
861 c$$$      iv(23)=0                                                          
862 c$$$c      iv(23)=1                  ! DEBUG
863 c$$$* 1 means to print initial x and d                                      
864 c$$$      iv(24)=0                                                          
865 c$$$c      iv(24)=1                  ! DEBUG
866 c$$$* min val for v(radfac) default is 0.1                                  
867 c$$$      v(24)=0.1D0                                                       
868 c$$$* max val for v(radfac) default is 4.0                                  
869 c$$$      v(25)=2.0D0                                                       
870 c$$$c     v(25)=4.0D0                                                       
871 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
872 c$$$* the sumsl default is 0.1                                              
873 c$$$      v(26)=0.1D0
874 c$$$* false conv if (act fnctn decrease) .lt. v(34)                         
875 c$$$* the sumsl default is 100*machep                                       
876 c$$$      v(34)=v(34)/100.0D0                                               
877 c$$$* absolute convergence                                                  
878 c$$$      if (tolf.eq.0.0D0) tolf=1.0D-4
879 c$$$      v(31)=tolf
880 c$$$* relative convergence                                                  
881 c$$$      if (rtolf.eq.0.0D0) rtolf=1.0D-1
882 c$$$      v(32)=rtolf
883 c$$$* controls initial step size                                            
884 c$$$       v(35)=1.0D-1                                                    
885 c$$$* large vals of d correspond to small components of step                
886 c$$$      do i=1,nphi
887 c$$$        d(i)=1.0D-1
888 c$$$      enddo
889 c$$$      do i=nphi+1,nvar
890 c$$$        d(i)=1.0D-1
891 c$$$      enddo
892 c$$$
893 c$$$      call geom_to_var(nvar,x)
894 c$$$      IF (mask_r) THEN
895 c$$$        do i=1,nres             ! Just in case...
896 c$$$          mask_phi(i)=0
897 c$$$          mask_theta(i)=0
898 c$$$        enddo
899 c$$$        call x2xx(x,xx,nvar_restr)
900 c$$$        call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
901 c$$$     &       iv,liv,lv,v,idum,rdum,fdum)      
902 c$$$        call xx2x(x,xx)
903 c$$$      ELSE
904 c$$$c     When minimizing ALL side-chains, etotal_sc is a little
905 c$$$c     faster if we don't set mask_r
906 c$$$        do i=1,nres
907 c$$$          mask_phi(i)=0
908 c$$$          mask_theta(i)=0
909 c$$$          mask_side(i)=1
910 c$$$        enddo
911 c$$$        call x2xx(x,xx,nvar_restr)
912 c$$$        call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
913 c$$$     &       iv,liv,lv,v,idum,rdum,fdum)      
914 c$$$        call xx2x(x,xx)
915 c$$$      ENDIF
916 c$$$      call var_to_geom(nvar,x)
917 c$$$      call chainbuild_sc
918 c$$$      etot_sc=v(10)                                                      
919 c$$$      iretcode=iv(1)
920 c$$$      nfun=iv(6)
921 c$$$      return  
922 c$$$      end  
923 c$$$
924 c$$$C--------------------------------------------------------------------------
925 c$$$
926 c$$$      subroutine chainbuild_sc
927 c$$$      implicit none
928 c$$$      include 'DIMENSIONS'
929 c$$$      include 'COMMON.VAR'
930 c$$$      include 'COMMON.INTERACT'
931 c$$$
932 c$$$c     Local variables
933 c$$$      integer i
934 c$$$
935 c$$$
936 c$$$      do i=nnt,nct
937 c$$$        if (.not.mask_r .or. mask_side(i).eq.1) then
938 c$$$          call locate_side_chain(i)
939 c$$$        endif
940 c$$$      enddo
941 c$$$
942 c$$$      return
943 c$$$      end
944 c$$$
945 c$$$C--------------------------------------------------------------------------
946 c$$$
947 c$$$      subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)  
948 c$$$      implicit none
949 c$$$
950 c$$$c     Includes
951 c$$$      include 'DIMENSIONS'
952 c$$$      include 'COMMON.DERIV'
953 c$$$      include 'COMMON.VAR'
954 c$$$      include 'COMMON.MINIM'
955 c$$$      include 'COMMON.IOUNITS'
956 c$$$
957 c$$$c     Input arguments
958 c$$$      integer n
959 c$$$      double precision x(maxvar)
960 c$$$      double precision ufparm
961 c$$$      external ufparm
962 c$$$
963 c$$$c     Input/Output arguments
964 c$$$      integer nf
965 c$$$      integer uiparm(1)
966 c$$$      double precision urparm(1)
967 c$$$
968 c$$$c     Output arguments
969 c$$$      double precision f
970 c$$$
971 c$$$c     Local variables
972 c$$$      double precision energia(0:n_ene)
973 c$$$#ifdef OSF
974 c$$$c     Variables used to intercept NaNs
975 c$$$      double precision x_sum
976 c$$$      integer i_NAN
977 c$$$#endif
978 c$$$
979 c$$$
980 c$$$      nfl=nf
981 c$$$      icg=mod(nf,2)+1
982 c$$$
983 c$$$#ifdef OSF
984 c$$$c     Intercept NaNs in the coordinates, before calling etotal_sc
985 c$$$      x_sum=0.D0
986 c$$$      do i_NAN=1,n
987 c$$$        x_sum=x_sum+x(i_NAN)
988 c$$$      enddo
989 c$$$c     Calculate the energy only if the coordinates are ok
990 c$$$      if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
991 c$$$        write(iout,*)"   *** func_restr_sc : Found NaN in coordinates"
992 c$$$        f=1.0D+77
993 c$$$        nf=0
994 c$$$      else
995 c$$$#endif
996 c$$$
997 c$$$      call var_to_geom_restr(n,x)
998 c$$$      call zerograd
999 c$$$      call chainbuild_sc
1000 c$$$      call etotal_sc(energia(0))
1001 c$$$      f=energia(0)
1002 c$$$      if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
1003 c$$$
1004 c$$$#ifdef OSF
1005 c$$$      endif
1006 c$$$#endif
1007 c$$$
1008 c$$$      return                                                            
1009 c$$$      end                                                               
1010 c$$$
1011 c$$$c-------------------------------------------------------
1012 c$$$
1013 c$$$      subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
1014 c$$$      implicit none
1015 c$$$
1016 c$$$c     Includes
1017 c$$$      include 'DIMENSIONS'
1018 c$$$      include 'COMMON.CHAIN'
1019 c$$$      include 'COMMON.DERIV'
1020 c$$$      include 'COMMON.VAR'
1021 c$$$      include 'COMMON.INTERACT'
1022 c$$$      include 'COMMON.MINIM'
1023 c$$$
1024 c$$$c     Input arguments
1025 c$$$      integer n
1026 c$$$      double precision x(maxvar)
1027 c$$$      double precision ufparm
1028 c$$$      external ufparm
1029 c$$$
1030 c$$$c     Input/Output arguments
1031 c$$$      integer nf
1032 c$$$      integer uiparm(1)
1033 c$$$      double precision urparm(1)
1034 c$$$
1035 c$$$c     Output arguments
1036 c$$$      double precision g(maxvar)
1037 c$$$
1038 c$$$c     Local variables
1039 c$$$      double precision f,gphii,gthetai,galphai,gomegai
1040 c$$$      integer ig,ind,i,j,k,igall,ij
1041 c$$$
1042 c$$$
1043 c$$$      icg=mod(nf,2)+1
1044 c$$$      if (nf-nfl+1) 20,30,40
1045 c$$$   20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1046 c$$$c     write (iout,*) 'grad 20'
1047 c$$$      if (nf.eq.0) return
1048 c$$$      goto 40
1049 c$$$   30 call var_to_geom_restr(n,x)
1050 c$$$      call chainbuild_sc
1051 c$$$C
1052 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1053 c$$$C
1054 c$$$   40 call cartder
1055 c$$$C
1056 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1057 c$$$C
1058 c$$$
1059 c$$$      ig=0
1060 c$$$      ind=nres-2
1061 c$$$      do i=2,nres-2
1062 c$$$       IF (mask_phi(i+2).eq.1) THEN
1063 c$$$        gphii=0.0D0
1064 c$$$        do j=i+1,nres-1
1065 c$$$          ind=ind+1
1066 c$$$          do k=1,3
1067 c$$$            gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1068 c$$$            gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1069 c$$$          enddo
1070 c$$$        enddo
1071 c$$$        ig=ig+1
1072 c$$$        g(ig)=gphii
1073 c$$$       ELSE
1074 c$$$        ind=ind+nres-1-i
1075 c$$$       ENDIF
1076 c$$$      enddo                                        
1077 c$$$
1078 c$$$
1079 c$$$      ind=0
1080 c$$$      do i=1,nres-2
1081 c$$$       IF (mask_theta(i+2).eq.1) THEN
1082 c$$$        ig=ig+1
1083 c$$$    gthetai=0.0D0
1084 c$$$    do j=i+1,nres-1
1085 c$$$          ind=ind+1
1086 c$$$      do k=1,3
1087 c$$$            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1088 c$$$            gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1089 c$$$          enddo
1090 c$$$        enddo
1091 c$$$        g(ig)=gthetai
1092 c$$$       ELSE
1093 c$$$        ind=ind+nres-1-i
1094 c$$$       ENDIF
1095 c$$$      enddo
1096 c$$$
1097 c$$$      do i=2,nres-1
1098 c$$$    if (itype(i).ne.10) then
1099 c$$$         IF (mask_side(i).eq.1) THEN
1100 c$$$          ig=ig+1
1101 c$$$          galphai=0.0D0
1102 c$$$      do k=1,3
1103 c$$$        galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1104 c$$$          enddo
1105 c$$$          g(ig)=galphai
1106 c$$$         ENDIF
1107 c$$$        endif
1108 c$$$      enddo
1109 c$$$
1110 c$$$      
1111 c$$$      do i=2,nres-1
1112 c$$$        if (itype(i).ne.10) then
1113 c$$$         IF (mask_side(i).eq.1) THEN
1114 c$$$          ig=ig+1
1115 c$$$      gomegai=0.0D0
1116 c$$$      do k=1,3
1117 c$$$        gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1118 c$$$          enddo
1119 c$$$      g(ig)=gomegai
1120 c$$$         ENDIF
1121 c$$$        endif
1122 c$$$      enddo
1123 c$$$
1124 c$$$C
1125 c$$$C Add the components corresponding to local energy terms.
1126 c$$$C
1127 c$$$
1128 c$$$      ig=0
1129 c$$$      igall=0
1130 c$$$      do i=4,nres
1131 c$$$        igall=igall+1
1132 c$$$        if (mask_phi(i).eq.1) then
1133 c$$$          ig=ig+1
1134 c$$$          g(ig)=g(ig)+gloc(igall,icg)
1135 c$$$        endif
1136 c$$$      enddo
1137 c$$$
1138 c$$$      do i=3,nres
1139 c$$$        igall=igall+1
1140 c$$$        if (mask_theta(i).eq.1) then
1141 c$$$          ig=ig+1
1142 c$$$          g(ig)=g(ig)+gloc(igall,icg)
1143 c$$$        endif
1144 c$$$      enddo
1145 c$$$     
1146 c$$$      do ij=1,2
1147 c$$$      do i=2,nres-1
1148 c$$$        if (itype(i).ne.10) then
1149 c$$$          igall=igall+1
1150 c$$$          if (mask_side(i).eq.1) then
1151 c$$$            ig=ig+1
1152 c$$$            g(ig)=g(ig)+gloc(igall,icg)
1153 c$$$          endif
1154 c$$$        endif
1155 c$$$      enddo
1156 c$$$      enddo
1157 c$$$
1158 c$$$cd      do i=1,ig
1159 c$$$cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1160 c$$$cd      enddo
1161 c$$$
1162 c$$$      return
1163 c$$$      end
1164 c$$$
1165 c$$$C-----------------------------------------------------------------------------
1166 c$$$
1167 c$$$      subroutine etotal_sc(energy_sc)
1168 c$$$      implicit none
1169 c$$$
1170 c$$$c     Includes
1171 c$$$      include 'DIMENSIONS'
1172 c$$$      include 'COMMON.VAR'
1173 c$$$      include 'COMMON.INTERACT'
1174 c$$$      include 'COMMON.DERIV'
1175 c$$$      include 'COMMON.FFIELD'
1176 c$$$
1177 c$$$c     Output arguments
1178 c$$$      double precision energy_sc(0:n_ene)
1179 c$$$
1180 c$$$c     Local variables
1181 c$$$      double precision evdw,escloc
1182 c$$$      integer i,j
1183 c$$$
1184 c$$$
1185 c$$$      do i=1,n_ene
1186 c$$$        energy_sc(i)=0.0D0
1187 c$$$      enddo
1188 c$$$
1189 c$$$      if (mask_r) then
1190 c$$$        call egb_sc(evdw)
1191 c$$$        call esc_sc(escloc)
1192 c$$$      else
1193 c$$$        call egb(evdw)
1194 c$$$        call esc(escloc)
1195 c$$$      endif
1196 c$$$
1197 c$$$      if (evdw.eq.1.0D20) then
1198 c$$$        energy_sc(0)=evdw
1199 c$$$      else
1200 c$$$        energy_sc(0)=wsc*evdw+wscloc*escloc
1201 c$$$      endif
1202 c$$$      energy_sc(1)=evdw
1203 c$$$      energy_sc(12)=escloc
1204 c$$$
1205 c$$$C
1206 c$$$C Sum up the components of the Cartesian gradient.
1207 c$$$C
1208 c$$$      do i=1,nct
1209 c$$$        do j=1,3
1210 c$$$          gradx(j,i,icg)=wsc*gvdwx(j,i)
1211 c$$$        enddo
1212 c$$$      enddo
1213 c$$$
1214 c$$$      return
1215 c$$$      end
1216 c$$$
1217 c$$$C-----------------------------------------------------------------------------
1218 c$$$
1219 c$$$      subroutine egb_sc(evdw)
1220 c$$$C
1221 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1222 c$$$C assuming the Gay-Berne potential of interaction.
1223 c$$$C
1224 c$$$      implicit real*8 (a-h,o-z)
1225 c$$$      include 'DIMENSIONS'
1226 c$$$      include 'COMMON.GEO'
1227 c$$$      include 'COMMON.VAR'
1228 c$$$      include 'COMMON.LOCAL'
1229 c$$$      include 'COMMON.CHAIN'
1230 c$$$      include 'COMMON.DERIV'
1231 c$$$      include 'COMMON.NAMES'
1232 c$$$      include 'COMMON.INTERACT'
1233 c$$$      include 'COMMON.IOUNITS'
1234 c$$$      include 'COMMON.CALC'
1235 c$$$      include 'COMMON.CONTROL'
1236 c$$$      logical lprn
1237 c$$$      evdw=0.0D0
1238 c$$$      energy_dec=.false.
1239 c$$$c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1240 c$$$      evdw=0.0D0
1241 c$$$      lprn=.false.
1242 c$$$c     if (icall.eq.0) lprn=.false.
1243 c$$$      ind=0
1244 c$$$      do i=iatsc_s,iatsc_e
1245 c$$$        itypi=itype(i)
1246 c$$$        itypi1=itype(i+1)
1247 c$$$        xi=c(1,nres+i)
1248 c$$$        yi=c(2,nres+i)
1249 c$$$        zi=c(3,nres+i)
1250 c$$$        dxi=dc_norm(1,nres+i)
1251 c$$$        dyi=dc_norm(2,nres+i)
1252 c$$$        dzi=dc_norm(3,nres+i)
1253 c$$$c        dsci_inv=dsc_inv(itypi)
1254 c$$$        dsci_inv=vbld_inv(i+nres)
1255 c$$$c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1256 c$$$c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1257 c$$$C
1258 c$$$C Calculate SC interaction energy.
1259 c$$$C
1260 c$$$        do iint=1,nint_gr(i)
1261 c$$$          do j=istart(i,iint),iend(i,iint)
1262 c$$$          IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1263 c$$$            ind=ind+1
1264 c$$$            itypj=itype(j)
1265 c$$$c            dscj_inv=dsc_inv(itypj)
1266 c$$$            dscj_inv=vbld_inv(j+nres)
1267 c$$$c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1268 c$$$c     &       1.0d0/vbld(j+nres)
1269 c$$$c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1270 c$$$            sig0ij=sigma(itypi,itypj)
1271 c$$$            chi1=chi(itypi,itypj)
1272 c$$$            chi2=chi(itypj,itypi)
1273 c$$$            chi12=chi1*chi2
1274 c$$$            chip1=chip(itypi)
1275 c$$$            chip2=chip(itypj)
1276 c$$$            chip12=chip1*chip2
1277 c$$$            alf1=alp(itypi)
1278 c$$$            alf2=alp(itypj)
1279 c$$$            alf12=0.5D0*(alf1+alf2)
1280 c$$$C For diagnostics only!!!
1281 c$$$c           chi1=0.0D0
1282 c$$$c           chi2=0.0D0
1283 c$$$c           chi12=0.0D0
1284 c$$$c           chip1=0.0D0
1285 c$$$c           chip2=0.0D0
1286 c$$$c           chip12=0.0D0
1287 c$$$c           alf1=0.0D0
1288 c$$$c           alf2=0.0D0
1289 c$$$c           alf12=0.0D0
1290 c$$$            xj=c(1,nres+j)-xi
1291 c$$$            yj=c(2,nres+j)-yi
1292 c$$$            zj=c(3,nres+j)-zi
1293 c$$$            dxj=dc_norm(1,nres+j)
1294 c$$$            dyj=dc_norm(2,nres+j)
1295 c$$$            dzj=dc_norm(3,nres+j)
1296 c$$$c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1297 c$$$c            write (iout,*) "j",j," dc_norm",
1298 c$$$c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1299 c$$$            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1300 c$$$            rij=dsqrt(rrij)
1301 c$$$C Calculate angle-dependent terms of energy and contributions to their
1302 c$$$C derivatives.
1303 c$$$            call sc_angular
1304 c$$$            sigsq=1.0D0/sigsq
1305 c$$$            sig=sig0ij*dsqrt(sigsq)
1306 c$$$            rij_shift=1.0D0/rij-sig+sig0ij
1307 c$$$c for diagnostics; uncomment
1308 c$$$c            rij_shift=1.2*sig0ij
1309 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1310 c$$$            if (rij_shift.le.0.0D0) then
1311 c$$$              evdw=1.0D20
1312 c$$$cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1313 c$$$cd     &        restyp(itypi),i,restyp(itypj),j,
1314 c$$$cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
1315 c$$$              return
1316 c$$$            endif
1317 c$$$            sigder=-sig*sigsq
1318 c$$$c---------------------------------------------------------------
1319 c$$$            rij_shift=1.0D0/rij_shift 
1320 c$$$            fac=rij_shift**expon
1321 c$$$            e1=fac*fac*aa(itypi,itypj)
1322 c$$$            e2=fac*bb(itypi,itypj)
1323 c$$$            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1324 c$$$            eps2der=evdwij*eps3rt
1325 c$$$            eps3der=evdwij*eps2rt
1326 c$$$c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1327 c$$$c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1328 c$$$            evdwij=evdwij*eps2rt*eps3rt
1329 c$$$            evdw=evdw+evdwij
1330 c$$$            if (lprn) then
1331 c$$$            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1332 c$$$            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1333 c$$$            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1334 c$$$     &        restyp(itypi),i,restyp(itypj),j,
1335 c$$$     &        epsi,sigm,chi1,chi2,chip1,chip2,
1336 c$$$     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1337 c$$$     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1338 c$$$     &        evdwij
1339 c$$$            endif
1340 c$$$
1341 c$$$            if (energy_dec) write (iout,'(a6,2i,0pf7.3)') 
1342 c$$$     &                        'evdw',i,j,evdwij
1343 c$$$
1344 c$$$C Calculate gradient components.
1345 c$$$            e1=e1*eps1*eps2rt**2*eps3rt**2
1346 c$$$            fac=-expon*(e1+evdwij)*rij_shift
1347 c$$$            sigder=fac*sigder
1348 c$$$            fac=rij*fac
1349 c$$$c            fac=0.0d0
1350 c$$$C Calculate the radial part of the gradient
1351 c$$$            gg(1)=xj*fac
1352 c$$$            gg(2)=yj*fac
1353 c$$$            gg(3)=zj*fac
1354 c$$$C Calculate angular part of the gradient.
1355 c$$$            call sc_grad
1356 c$$$          ENDIF
1357 c$$$          enddo      ! j
1358 c$$$        enddo        ! iint
1359 c$$$      enddo          ! i
1360 c$$$      energy_dec=.false.
1361 c$$$      return
1362 c$$$      end
1363 c$$$
1364 c$$$c-----------------------------------------------------------------------------
1365 c$$$
1366 c$$$      subroutine esc_sc(escloc)
1367 c$$$C Calculate the local energy of a side chain and its derivatives in the
1368 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles 
1369 c$$$C ALPHA and OMEGA.
1370 c$$$      implicit real*8 (a-h,o-z)
1371 c$$$      include 'DIMENSIONS'
1372 c$$$      include 'COMMON.GEO'
1373 c$$$      include 'COMMON.LOCAL'
1374 c$$$      include 'COMMON.VAR'
1375 c$$$      include 'COMMON.INTERACT'
1376 c$$$      include 'COMMON.DERIV'
1377 c$$$      include 'COMMON.CHAIN'
1378 c$$$      include 'COMMON.IOUNITS'
1379 c$$$      include 'COMMON.NAMES'
1380 c$$$      include 'COMMON.FFIELD'
1381 c$$$      include 'COMMON.CONTROL'
1382 c$$$      double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1383 c$$$     &     ddersc0(3),ddummy(3),xtemp(3),temp(3)
1384 c$$$      common /sccalc/ time11,time12,time112,theti,it,nlobit
1385 c$$$      delta=0.02d0*pi
1386 c$$$      escloc=0.0D0
1387 c$$$c     write (iout,'(a)') 'ESC'
1388 c$$$      do i=loc_start,loc_end
1389 c$$$      IF (mask_side(i).eq.1) THEN
1390 c$$$        it=itype(i)
1391 c$$$        if (it.eq.10) goto 1
1392 c$$$        nlobit=nlob(it)
1393 c$$$c       print *,'i=',i,' it=',it,' nlobit=',nlobit
1394 c$$$c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1395 c$$$        theti=theta(i+1)-pipol
1396 c$$$        x(1)=dtan(theti)
1397 c$$$        x(2)=alph(i)
1398 c$$$        x(3)=omeg(i)
1399 c$$$
1400 c$$$        if (x(2).gt.pi-delta) then
1401 c$$$          xtemp(1)=x(1)
1402 c$$$          xtemp(2)=pi-delta
1403 c$$$          xtemp(3)=x(3)
1404 c$$$          call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1405 c$$$          xtemp(2)=pi
1406 c$$$          call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1407 c$$$          call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1408 c$$$     &        escloci,dersc(2))
1409 c$$$          call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1410 c$$$     &        ddersc0(1),dersc(1))
1411 c$$$          call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1412 c$$$     &        ddersc0(3),dersc(3))
1413 c$$$          xtemp(2)=pi-delta
1414 c$$$          call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1415 c$$$          xtemp(2)=pi
1416 c$$$          call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1417 c$$$          call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1418 c$$$     &            dersc0(2),esclocbi,dersc02)
1419 c$$$          call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1420 c$$$     &            dersc12,dersc01)
1421 c$$$          call splinthet(x(2),0.5d0*delta,ss,ssd)
1422 c$$$          dersc0(1)=dersc01
1423 c$$$          dersc0(2)=dersc02
1424 c$$$          dersc0(3)=0.0d0
1425 c$$$          do k=1,3
1426 c$$$            dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1427 c$$$          enddo
1428 c$$$          dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1429 c$$$c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1430 c$$$c    &             esclocbi,ss,ssd
1431 c$$$          escloci=ss*escloci+(1.0d0-ss)*esclocbi
1432 c$$$c         escloci=esclocbi
1433 c$$$c         write (iout,*) escloci
1434 c$$$        else if (x(2).lt.delta) then
1435 c$$$          xtemp(1)=x(1)
1436 c$$$          xtemp(2)=delta
1437 c$$$          xtemp(3)=x(3)
1438 c$$$          call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1439 c$$$          xtemp(2)=0.0d0
1440 c$$$          call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1441 c$$$          call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1442 c$$$     &        escloci,dersc(2))
1443 c$$$          call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1444 c$$$     &        ddersc0(1),dersc(1))
1445 c$$$          call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1446 c$$$     &        ddersc0(3),dersc(3))
1447 c$$$          xtemp(2)=delta
1448 c$$$          call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1449 c$$$          xtemp(2)=0.0d0
1450 c$$$          call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1451 c$$$          call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1452 c$$$     &            dersc0(2),esclocbi,dersc02)
1453 c$$$          call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1454 c$$$     &            dersc12,dersc01)
1455 c$$$          dersc0(1)=dersc01
1456 c$$$          dersc0(2)=dersc02
1457 c$$$          dersc0(3)=0.0d0
1458 c$$$          call splinthet(x(2),0.5d0*delta,ss,ssd)
1459 c$$$          do k=1,3
1460 c$$$            dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1461 c$$$          enddo
1462 c$$$          dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1463 c$$$c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1464 c$$$c    &             esclocbi,ss,ssd
1465 c$$$          escloci=ss*escloci+(1.0d0-ss)*esclocbi
1466 c$$$c         write (iout,*) escloci
1467 c$$$        else
1468 c$$$          call enesc(x,escloci,dersc,ddummy,.false.)
1469 c$$$        endif
1470 c$$$
1471 c$$$        escloc=escloc+escloci
1472 c$$$        if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1473 c$$$     &     'escloc',i,escloci
1474 c$$$c       write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1475 c$$$
1476 c$$$        gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1477 c$$$     &   wscloc*dersc(1)
1478 c$$$        gloc(ialph(i,1),icg)=wscloc*dersc(2)
1479 c$$$        gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1480 c$$$    1   continue
1481 c$$$      ENDIF
1482 c$$$      enddo
1483 c$$$      return
1484 c$$$      end
1485 c$$$
1486 c$$$C-----------------------------------------------------------------------------
1487 c$$$
1488 c$$$      subroutine egb_ij(i_sc,j_sc,evdw)
1489 c$$$C
1490 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1491 c$$$C assuming the Gay-Berne potential of interaction.
1492 c$$$C
1493 c$$$      implicit real*8 (a-h,o-z)
1494 c$$$      include 'DIMENSIONS'
1495 c$$$      include 'COMMON.GEO'
1496 c$$$      include 'COMMON.VAR'
1497 c$$$      include 'COMMON.LOCAL'
1498 c$$$      include 'COMMON.CHAIN'
1499 c$$$      include 'COMMON.DERIV'
1500 c$$$      include 'COMMON.NAMES'
1501 c$$$      include 'COMMON.INTERACT'
1502 c$$$      include 'COMMON.IOUNITS'
1503 c$$$      include 'COMMON.CALC'
1504 c$$$      include 'COMMON.CONTROL'
1505 c$$$      logical lprn
1506 c$$$      evdw=0.0D0
1507 c$$$      energy_dec=.false.
1508 c$$$c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1509 c$$$      evdw=0.0D0
1510 c$$$      lprn=.false.
1511 c$$$      ind=0
1512 c$$$c$$$      do i=iatsc_s,iatsc_e
1513 c$$$      i=i_sc
1514 c$$$        itypi=itype(i)
1515 c$$$        itypi1=itype(i+1)
1516 c$$$        xi=c(1,nres+i)
1517 c$$$        yi=c(2,nres+i)
1518 c$$$        zi=c(3,nres+i)
1519 c$$$        dxi=dc_norm(1,nres+i)
1520 c$$$        dyi=dc_norm(2,nres+i)
1521 c$$$        dzi=dc_norm(3,nres+i)
1522 c$$$c        dsci_inv=dsc_inv(itypi)
1523 c$$$        dsci_inv=vbld_inv(i+nres)
1524 c$$$c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1525 c$$$c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1526 c$$$C
1527 c$$$C Calculate SC interaction energy.
1528 c$$$C
1529 c$$$c$$$        do iint=1,nint_gr(i)
1530 c$$$c$$$          do j=istart(i,iint),iend(i,iint)
1531 c$$$        j=j_sc
1532 c$$$            ind=ind+1
1533 c$$$            itypj=itype(j)
1534 c$$$c            dscj_inv=dsc_inv(itypj)
1535 c$$$            dscj_inv=vbld_inv(j+nres)
1536 c$$$c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1537 c$$$c     &       1.0d0/vbld(j+nres)
1538 c$$$c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1539 c$$$            sig0ij=sigma(itypi,itypj)
1540 c$$$            chi1=chi(itypi,itypj)
1541 c$$$            chi2=chi(itypj,itypi)
1542 c$$$            chi12=chi1*chi2
1543 c$$$            chip1=chip(itypi)
1544 c$$$            chip2=chip(itypj)
1545 c$$$            chip12=chip1*chip2
1546 c$$$            alf1=alp(itypi)
1547 c$$$            alf2=alp(itypj)
1548 c$$$            alf12=0.5D0*(alf1+alf2)
1549 c$$$C For diagnostics only!!!
1550 c$$$c           chi1=0.0D0
1551 c$$$c           chi2=0.0D0
1552 c$$$c           chi12=0.0D0
1553 c$$$c           chip1=0.0D0
1554 c$$$c           chip2=0.0D0
1555 c$$$c           chip12=0.0D0
1556 c$$$c           alf1=0.0D0
1557 c$$$c           alf2=0.0D0
1558 c$$$c           alf12=0.0D0
1559 c$$$            xj=c(1,nres+j)-xi
1560 c$$$            yj=c(2,nres+j)-yi
1561 c$$$            zj=c(3,nres+j)-zi
1562 c$$$            dxj=dc_norm(1,nres+j)
1563 c$$$            dyj=dc_norm(2,nres+j)
1564 c$$$            dzj=dc_norm(3,nres+j)
1565 c$$$c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1566 c$$$c            write (iout,*) "j",j," dc_norm",
1567 c$$$c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1568 c$$$            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1569 c$$$            rij=dsqrt(rrij)
1570 c$$$C Calculate angle-dependent terms of energy and contributions to their
1571 c$$$C derivatives.
1572 c$$$            call sc_angular
1573 c$$$            sigsq=1.0D0/sigsq
1574 c$$$            sig=sig0ij*dsqrt(sigsq)
1575 c$$$            rij_shift=1.0D0/rij-sig+sig0ij
1576 c$$$c for diagnostics; uncomment
1577 c$$$c            rij_shift=1.2*sig0ij
1578 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1579 c$$$            if (rij_shift.le.0.0D0) then
1580 c$$$              evdw=1.0D20
1581 c$$$cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1582 c$$$cd     &        restyp(itypi),i,restyp(itypj),j,
1583 c$$$cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
1584 c$$$              return
1585 c$$$            endif
1586 c$$$            sigder=-sig*sigsq
1587 c$$$c---------------------------------------------------------------
1588 c$$$            rij_shift=1.0D0/rij_shift 
1589 c$$$            fac=rij_shift**expon
1590 c$$$            e1=fac*fac*aa(itypi,itypj)
1591 c$$$            e2=fac*bb(itypi,itypj)
1592 c$$$            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1593 c$$$            eps2der=evdwij*eps3rt
1594 c$$$            eps3der=evdwij*eps2rt
1595 c$$$c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1596 c$$$c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1597 c$$$            evdwij=evdwij*eps2rt*eps3rt
1598 c$$$            evdw=evdw+evdwij
1599 c$$$            if (lprn) then
1600 c$$$            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1601 c$$$            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1602 c$$$            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1603 c$$$     &        restyp(itypi),i,restyp(itypj),j,
1604 c$$$     &        epsi,sigm,chi1,chi2,chip1,chip2,
1605 c$$$     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1606 c$$$     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1607 c$$$     &        evdwij
1608 c$$$            endif
1609 c$$$
1610 c$$$            if (energy_dec) write (iout,'(a6,2i,0pf7.3)') 
1611 c$$$     &                        'evdw',i,j,evdwij
1612 c$$$
1613 c$$$C Calculate gradient components.
1614 c$$$            e1=e1*eps1*eps2rt**2*eps3rt**2
1615 c$$$            fac=-expon*(e1+evdwij)*rij_shift
1616 c$$$            sigder=fac*sigder
1617 c$$$            fac=rij*fac
1618 c$$$c            fac=0.0d0
1619 c$$$C Calculate the radial part of the gradient
1620 c$$$            gg(1)=xj*fac
1621 c$$$            gg(2)=yj*fac
1622 c$$$            gg(3)=zj*fac
1623 c$$$C Calculate angular part of the gradient.
1624 c$$$            call sc_grad
1625 c$$$c$$$          enddo      ! j
1626 c$$$c$$$        enddo        ! iint
1627 c$$$c$$$      enddo          ! i
1628 c$$$      energy_dec=.false.
1629 c$$$      return
1630 c$$$      end
1631 c$$$
1632 c$$$C-----------------------------------------------------------------------------
1633 c$$$
1634 c$$$      subroutine perturb_side_chain(i,angle)
1635 c$$$      implicit none
1636 c$$$
1637 c$$$c     Includes
1638 c$$$      include 'DIMENSIONS'
1639 c$$$      include 'COMMON.CHAIN'
1640 c$$$      include 'COMMON.GEO'
1641 c$$$      include 'COMMON.VAR'
1642 c$$$      include 'COMMON.LOCAL'
1643 c$$$      include 'COMMON.IOUNITS'
1644 c$$$
1645 c$$$c     External functions
1646 c$$$      external ran_number
1647 c$$$      double precision ran_number
1648 c$$$
1649 c$$$c     Input arguments
1650 c$$$      integer i
1651 c$$$      double precision angle    ! In degrees
1652 c$$$
1653 c$$$c     Local variables
1654 c$$$      integer i_sc
1655 c$$$      double precision rad_ang,rand_v(3),length,cost,sint
1656 c$$$
1657 c$$$
1658 c$$$      i_sc=i+nres
1659 c$$$      rad_ang=angle*deg2rad
1660 c$$$
1661 c$$$      length=0.0
1662 c$$$      do while (length.lt.0.01)
1663 c$$$        rand_v(1)=ran_number(0.01D0,1.0D0)
1664 c$$$        rand_v(2)=ran_number(0.01D0,1.0D0)
1665 c$$$        rand_v(3)=ran_number(0.01D0,1.0D0)
1666 c$$$        length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1667 c$$$     +       rand_v(3)*rand_v(3)
1668 c$$$        length=sqrt(length)
1669 c$$$        rand_v(1)=rand_v(1)/length
1670 c$$$        rand_v(2)=rand_v(2)/length
1671 c$$$        rand_v(3)=rand_v(3)/length
1672 c$$$        cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1673 c$$$     +       rand_v(3)*dc_norm(3,i_sc)
1674 c$$$        length=1.0D0-cost*cost
1675 c$$$        if (length.lt.0.0D0) length=0.0D0
1676 c$$$        length=sqrt(length)
1677 c$$$        rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1678 c$$$        rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1679 c$$$        rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1680 c$$$      enddo
1681 c$$$      rand_v(1)=rand_v(1)/length
1682 c$$$      rand_v(2)=rand_v(2)/length
1683 c$$$      rand_v(3)=rand_v(3)/length
1684 c$$$
1685 c$$$      cost=dcos(rad_ang)
1686 c$$$      sint=dsin(rad_ang)
1687 c$$$      dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1688 c$$$      dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1689 c$$$      dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1690 c$$$      dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1691 c$$$      dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1692 c$$$      dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1693 c$$$      c(1,i_sc)=c(1,i)+dc(1,i_sc)
1694 c$$$      c(2,i_sc)=c(2,i)+dc(2,i_sc)
1695 c$$$      c(3,i_sc)=c(3,i)+dc(3,i_sc)
1696 c$$$
1697 c$$$      call chainbuild_cart
1698 c$$$
1699 c$$$      return
1700 c$$$      end
1701 c$$$
1702 c$$$c----------------------------------------------------------------------------
1703 c$$$
1704 c$$$      subroutine ss_relax3(i_in,j_in)
1705 c$$$      implicit none
1706 c$$$
1707 c$$$c     Includes
1708 c$$$      include 'DIMENSIONS'
1709 c$$$      include 'COMMON.VAR'
1710 c$$$      include 'COMMON.CHAIN'
1711 c$$$      include 'COMMON.IOUNITS'
1712 c$$$      include 'COMMON.INTERACT'
1713 c$$$
1714 c$$$c     External functions
1715 c$$$      external ran_number
1716 c$$$      double precision ran_number
1717 c$$$
1718 c$$$c     Input arguments
1719 c$$$      integer i_in,j_in
1720 c$$$
1721 c$$$c     Local variables
1722 c$$$      double precision energy_sc(0:n_ene),etot
1723 c$$$      double precision org_dc(3),org_dc_norm(3),org_c(3)
1724 c$$$      double precision ang_pert,rand_fact,exp_fact,beta
1725 c$$$      integer n,i_pert,i
1726 c$$$      logical notdone
1727 c$$$
1728 c$$$
1729 c$$$      beta=1.0D0
1730 c$$$
1731 c$$$      mask_r=.true.
1732 c$$$      do i=nnt,nct
1733 c$$$        mask_side(i)=0
1734 c$$$      enddo
1735 c$$$      mask_side(i_in)=1
1736 c$$$      mask_side(j_in)=1
1737 c$$$
1738 c$$$      call etotal_sc(energy_sc)
1739 c$$$      etot=energy_sc(0)
1740 c$$$c      write(iout,'(a,3d15.5)')"     SS_MC_START ",energy_sc(0),
1741 c$$$c     +     energy_sc(1),energy_sc(12)
1742 c$$$
1743 c$$$      notdone=.true.
1744 c$$$      n=0
1745 c$$$      do while (notdone)
1746 c$$$        if (mod(n,2).eq.0) then
1747 c$$$          i_pert=i_in
1748 c$$$        else
1749 c$$$          i_pert=j_in
1750 c$$$        endif
1751 c$$$        n=n+1
1752 c$$$
1753 c$$$        do i=1,3
1754 c$$$          org_dc(i)=dc(i,i_pert+nres)
1755 c$$$          org_dc_norm(i)=dc_norm(i,i_pert+nres)
1756 c$$$          org_c(i)=c(i,i_pert+nres)
1757 c$$$        enddo
1758 c$$$        ang_pert=ran_number(0.0D0,3.0D0)
1759 c$$$        call perturb_side_chain(i_pert,ang_pert)
1760 c$$$        call etotal_sc(energy_sc)
1761 c$$$        exp_fact=exp(beta*(etot-energy_sc(0)))
1762 c$$$        rand_fact=ran_number(0.0D0,1.0D0)
1763 c$$$        if (rand_fact.lt.exp_fact) then
1764 c$$$c          write(iout,'(a,3d15.5)')"     SS_MC_ACCEPT ",energy_sc(0),
1765 c$$$c     +     energy_sc(1),energy_sc(12)
1766 c$$$          etot=energy_sc(0)
1767 c$$$        else
1768 c$$$c          write(iout,'(a,3d15.5)')"     SS_MC_REJECT ",energy_sc(0),
1769 c$$$c     +     energy_sc(1),energy_sc(12)
1770 c$$$          do i=1,3
1771 c$$$            dc(i,i_pert+nres)=org_dc(i)
1772 c$$$            dc_norm(i,i_pert+nres)=org_dc_norm(i)
1773 c$$$            c(i,i_pert+nres)=org_c(i)
1774 c$$$          enddo
1775 c$$$        endif
1776 c$$$
1777 c$$$        if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1778 c$$$      enddo
1779 c$$$
1780 c$$$      mask_r=.false.
1781 c$$$
1782 c$$$      return
1783 c$$$      end
1784 c$$$
1785 c$$$c----------------------------------------------------------------------------
1786 c$$$
1787 c$$$      subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1788 c$$$      implicit none
1789 c$$$      include 'DIMENSIONS'
1790 c$$$      integer liv,lv
1791 c$$$      parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2)) 
1792 c$$$*********************************************************************
1793 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1794 c$$$* the calling subprogram.                                           *     
1795 c$$$* when d(i)=1.0, then v(35) is the length of the initial step,      *     
1796 c$$$* calculated in the usual pythagorean way.                          *     
1797 c$$$* absolute convergence occurs when the function is within v(31) of  *     
1798 c$$$* zero. unless you know the minimum value in advance, abs convg     *     
1799 c$$$* is probably not useful.                                           *     
1800 c$$$* relative convergence is when the model predicts that the function *   
1801 c$$$* will decrease by less than v(32)*abs(fun).                        *   
1802 c$$$*********************************************************************
1803 c$$$      include 'COMMON.IOUNITS'
1804 c$$$      include 'COMMON.VAR'
1805 c$$$      include 'COMMON.GEO'
1806 c$$$      include 'COMMON.MINIM'
1807 c$$$      include 'COMMON.CHAIN'
1808 c$$$
1809 c$$$      double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1810 c$$$      common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1811 c$$$     +     orig_ss_dist(maxres2,maxres2)
1812 c$$$
1813 c$$$      double precision etot
1814 c$$$      integer iretcode,nfun,i_in,j_in
1815 c$$$
1816 c$$$      external dist
1817 c$$$      double precision dist
1818 c$$$      external ss_func,fdum
1819 c$$$      double precision ss_func,fdum
1820 c$$$
1821 c$$$      integer iv(liv),uiparm(2)
1822 c$$$      double precision v(lv),x(maxres6),d(maxres6),rdum
1823 c$$$      integer i,j,k
1824 c$$$
1825 c$$$
1826 c$$$      call deflt(2,iv,liv,lv,v)                                         
1827 c$$$* 12 means fresh start, dont call deflt                                 
1828 c$$$      iv(1)=12                                                          
1829 c$$$* max num of fun calls                                                  
1830 c$$$      if (maxfun.eq.0) maxfun=500
1831 c$$$      iv(17)=maxfun
1832 c$$$* max num of iterations                                                 
1833 c$$$      if (maxmin.eq.0) maxmin=1000
1834 c$$$      iv(18)=maxmin
1835 c$$$* controls output                                                       
1836 c$$$      iv(19)=2                                                          
1837 c$$$* selects output unit                                                   
1838 c$$$c      iv(21)=iout                                                       
1839 c$$$      iv(21)=0
1840 c$$$* 1 means to print out result                                           
1841 c$$$      iv(22)=0                                                          
1842 c$$$* 1 means to print out summary stats                                    
1843 c$$$      iv(23)=0                                                          
1844 c$$$* 1 means to print initial x and d                                      
1845 c$$$      iv(24)=0                                                          
1846 c$$$* min val for v(radfac) default is 0.1                                  
1847 c$$$      v(24)=0.1D0                                                       
1848 c$$$* max val for v(radfac) default is 4.0                                  
1849 c$$$      v(25)=2.0D0                                                       
1850 c$$$c     v(25)=4.0D0                                                       
1851 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
1852 c$$$* the sumsl default is 0.1                                              
1853 c$$$      v(26)=0.1D0
1854 c$$$* false conv if (act fnctn decrease) .lt. v(34)                         
1855 c$$$* the sumsl default is 100*machep                                       
1856 c$$$      v(34)=v(34)/100.0D0                                               
1857 c$$$* absolute convergence                                                  
1858 c$$$      if (tolf.eq.0.0D0) tolf=1.0D-4
1859 c$$$      v(31)=tolf
1860 c$$$      v(31)=1.0D-1
1861 c$$$* relative convergence                                                  
1862 c$$$      if (rtolf.eq.0.0D0) rtolf=1.0D-4
1863 c$$$      v(32)=rtolf
1864 c$$$      v(32)=1.0D-1
1865 c$$$* controls initial step size                                            
1866 c$$$      v(35)=1.0D-1
1867 c$$$* large vals of d correspond to small components of step                
1868 c$$$      do i=1,6*nres
1869 c$$$        d(i)=1.0D0
1870 c$$$      enddo
1871 c$$$
1872 c$$$      do i=0,2*nres
1873 c$$$        do j=1,3
1874 c$$$          orig_ss_dc(j,i)=dc(j,i)
1875 c$$$        enddo
1876 c$$$      enddo
1877 c$$$      call geom_to_var(nvar,orig_ss_var)
1878 c$$$
1879 c$$$      do i=1,nres
1880 c$$$        do j=i,nres
1881 c$$$          orig_ss_dist(j,i)=dist(j,i)
1882 c$$$          orig_ss_dist(j+nres,i)=dist(j+nres,i)
1883 c$$$          orig_ss_dist(j,i+nres)=dist(j,i+nres)
1884 c$$$          orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1885 c$$$        enddo
1886 c$$$      enddo
1887 c$$$
1888 c$$$      k=0
1889 c$$$      do i=1,nres-1
1890 c$$$        do j=1,3
1891 c$$$          k=k+1
1892 c$$$          x(k)=dc(j,i)
1893 c$$$        enddo
1894 c$$$      enddo
1895 c$$$      do i=2,nres-1
1896 c$$$        if (ialph(i,1).gt.0) then
1897 c$$$        do j=1,3
1898 c$$$          k=k+1
1899 c$$$          x(k)=dc(j,i+nres)
1900 c$$$        enddo
1901 c$$$        endif
1902 c$$$      enddo
1903 c$$$
1904 c$$$      uiparm(1)=i_in
1905 c$$$      uiparm(2)=j_in
1906 c$$$      call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1907 c$$$      etot=v(10)
1908 c$$$      iretcode=iv(1)
1909 c$$$      nfun=iv(6)+iv(30)
1910 c$$$
1911 c$$$      k=0
1912 c$$$      do i=1,nres-1
1913 c$$$        do j=1,3
1914 c$$$          k=k+1
1915 c$$$          dc(j,i)=x(k)
1916 c$$$        enddo
1917 c$$$      enddo
1918 c$$$      do i=2,nres-1
1919 c$$$        if (ialph(i,1).gt.0) then
1920 c$$$        do j=1,3
1921 c$$$          k=k+1
1922 c$$$          dc(j,i+nres)=x(k)
1923 c$$$        enddo
1924 c$$$        endif
1925 c$$$      enddo
1926 c$$$      call chainbuild_cart
1927 c$$$
1928 c$$$      return  
1929 c$$$      end  
1930 c$$$
1931 c$$$C-----------------------------------------------------------------------------
1932 c$$$
1933 c$$$      subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)  
1934 c$$$      implicit none
1935 c$$$      include 'DIMENSIONS'
1936 c$$$      include 'COMMON.DERIV'
1937 c$$$      include 'COMMON.IOUNITS'
1938 c$$$      include 'COMMON.VAR'
1939 c$$$      include 'COMMON.CHAIN'
1940 c$$$      include 'COMMON.INTERACT'
1941 c$$$      include 'COMMON.SBRIDGE'
1942 c$$$
1943 c$$$      double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1944 c$$$      common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1945 c$$$     +     orig_ss_dist(maxres2,maxres2)
1946 c$$$
1947 c$$$      integer n
1948 c$$$      double precision x(maxres6)
1949 c$$$      integer nf
1950 c$$$      double precision f
1951 c$$$      integer uiparm(2)
1952 c$$$      real*8 urparm(1)
1953 c$$$      external ufparm
1954 c$$$      double precision ufparm
1955 c$$$
1956 c$$$      external dist
1957 c$$$      double precision dist
1958 c$$$
1959 c$$$      integer i,j,k,ss_i,ss_j
1960 c$$$      double precision tempf,var(maxvar)
1961 c$$$
1962 c$$$
1963 c$$$      ss_i=uiparm(1)
1964 c$$$      ss_j=uiparm(2)
1965 c$$$      f=0.0D0
1966 c$$$
1967 c$$$      k=0
1968 c$$$      do i=1,nres-1
1969 c$$$        do j=1,3
1970 c$$$          k=k+1
1971 c$$$          dc(j,i)=x(k)
1972 c$$$        enddo
1973 c$$$      enddo
1974 c$$$      do i=2,nres-1
1975 c$$$        if (ialph(i,1).gt.0) then
1976 c$$$        do j=1,3
1977 c$$$          k=k+1
1978 c$$$          dc(j,i+nres)=x(k)
1979 c$$$        enddo
1980 c$$$        endif
1981 c$$$      enddo
1982 c$$$      call chainbuild_cart
1983 c$$$
1984 c$$$      call geom_to_var(nvar,var)
1985 c$$$
1986 c$$$c     Constraints on all angles
1987 c$$$      do i=1,nvar
1988 c$$$        tempf=var(i)-orig_ss_var(i)
1989 c$$$        f=f+tempf*tempf
1990 c$$$      enddo
1991 c$$$
1992 c$$$c     Constraints on all distances
1993 c$$$      do i=1,nres-1
1994 c$$$        if (i.gt.1) then
1995 c$$$          tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
1996 c$$$          f=f+tempf*tempf
1997 c$$$        endif
1998 c$$$        do j=i+1,nres
1999 c$$$          tempf=dist(j,i)-orig_ss_dist(j,i)
2000 c$$$          if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
2001 c$$$          tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
2002 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
2003 c$$$          tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
2004 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
2005 c$$$          tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
2006 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
2007 c$$$        enddo
2008 c$$$      enddo
2009 c$$$
2010 c$$$c     Constraints for the relevant CYS-CYS
2011 c$$$      tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
2012 c$$$      f=f+tempf*tempf
2013 c$$$CCCCCCCCCCCCCCCCC      ADD SOME ANGULAR STUFF
2014 c$$$
2015 c$$$c$$$      if (nf.ne.nfl) then
2016 c$$$c$$$        write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
2017 c$$$c$$$     +       f,dist(5+nres,14+nres)
2018 c$$$c$$$      endif
2019 c$$$
2020 c$$$      nfl=nf
2021 c$$$
2022 c$$$      return                                                            
2023 c$$$      end                                                               
2024 c$$$
2025 c$$$C-----------------------------------------------------------------------------
2026 c$$$C-----------------------------------------------------------------------------
2027          subroutine triple_ssbond_ene(resi,resj,resk,eij)
2028       include 'DIMENSIONS'
2029       include 'COMMON.SBRIDGE'
2030       include 'COMMON.CHAIN'
2031       include 'COMMON.DERIV'
2032       include 'COMMON.LOCAL'
2033       include 'COMMON.INTERACT'
2034       include 'COMMON.VAR'
2035       include 'COMMON.IOUNITS'
2036       include 'COMMON.CALC'
2037 #ifndef CLUST
2038 #ifndef WHAM
2039 C      include 'COMMON.MD'
2040 #endif
2041 #endif
2042
2043 c     External functions
2044       double precision h_base
2045       external h_base
2046
2047 c     Input arguments
2048       integer resi,resj,resk
2049
2050 c     Output arguments
2051       double precision eij,eij1,eij2,eij3
2052
2053 c     Local variables
2054       logical havebond
2055 c      integer itypi,itypj,k,l
2056       double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
2057       double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
2058       double precision xik,yik,zik,xjk,yjk,zjk
2059       double precision sig0ij,ljd,sig,fac,e1,e2
2060       double precision dcosom1(3),dcosom2(3),ed
2061       double precision pom1,pom2
2062       double precision ljA,ljB,ljXs
2063       double precision d_ljB(1:3)
2064       double precision ssA,ssB,ssC,ssXs
2065       double precision ssxm,ljxm,ssm,ljm
2066       double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
2067
2068       i=resi
2069       j=resj
2070       k=resk
2071 C      write(iout,*) resi,resj,resk
2072       itypi=itype(i)
2073       dxi=dc_norm(1,nres+i)
2074       dyi=dc_norm(2,nres+i)
2075       dzi=dc_norm(3,nres+i)
2076       dsci_inv=vbld_inv(i+nres)
2077       xi=c(1,nres+i)
2078       yi=c(2,nres+i)
2079       zi=c(3,nres+i)
2080
2081       itypj=itype(j)
2082       xj=c(1,nres+j)
2083       yj=c(2,nres+j)
2084       zj=c(3,nres+j)
2085       
2086       dxj=dc_norm(1,nres+j)
2087       dyj=dc_norm(2,nres+j)
2088       dzj=dc_norm(3,nres+j)
2089       dscj_inv=vbld_inv(j+nres)
2090       itypk=itype(k)
2091       xk=c(1,nres+k)
2092       yk=c(2,nres+k)
2093       zk=c(3,nres+k)
2094       
2095       dxk=dc_norm(1,nres+k)
2096       dyk=dc_norm(2,nres+k)
2097       dzk=dc_norm(3,nres+k)
2098       dscj_inv=vbld_inv(k+nres)
2099       xij=xj-xi
2100       xik=xk-xi
2101       xjk=xk-xj
2102       yij=yj-yi
2103       yik=yk-yi
2104       yjk=yk-yj
2105       zij=zj-zi
2106       zik=zk-zi
2107       zjk=zk-zj
2108       rrij=(xij*xij+yij*yij+zij*zij)
2109       rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
2110       rrik=(xik*xik+yik*yik+zik*zik)
2111       rik=dsqrt(rrik)
2112       rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
2113       rjk=dsqrt(rrjk)
2114 C there are three combination of distances for each trisulfide bonds
2115 C The first case the ith atom is the center
2116 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
2117 C distance y is second distance the a,b,c,d are parameters derived for
2118 C this problem d parameter was set as a penalty currenlty set to 1.
2119       eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
2120 C second case jth atom is center
2121       eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
2122 C the third case kth atom is the center
2123       eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
2124 C      eij2=0.0
2125 C      eij3=0.0
2126 C      eij1=0.0
2127       eij=eij1+eij2+eij3
2128 C      write(iout,*)i,j,k,eij
2129 C The energy penalty calculated now time for the gradient part 
2130 C derivative over rij
2131       fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2132      &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))  
2133             gg(1)=xij*fac/rij
2134             gg(2)=yij*fac/rij
2135             gg(3)=zij*fac/rij
2136       do m=1,3
2137         gvdwx(m,i)=gvdwx(m,i)-gg(m)
2138         gvdwx(m,j)=gvdwx(m,j)+gg(m)
2139       enddo
2140       do l=1,3
2141         gvdwc(l,i)=gvdwc(l,i)-gg(l)
2142         gvdwc(l,j)=gvdwc(l,j)+gg(l)
2143       enddo
2144 C now derivative over rik
2145       fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
2146      &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2147             gg(1)=xik*fac/rik
2148             gg(2)=yik*fac/rik
2149             gg(3)=zik*fac/rik
2150       do m=1,3
2151         gvdwx(m,i)=gvdwx(m,i)-gg(m)
2152         gvdwx(m,k)=gvdwx(m,k)+gg(m)
2153       enddo
2154       do l=1,3
2155         gvdwc(l,i)=gvdwc(l,i)-gg(l)
2156         gvdwc(l,k)=gvdwc(l,k)+gg(l)
2157       enddo
2158 C now derivative over rjk
2159       fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
2160      &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
2161             gg(1)=xjk*fac/rjk
2162             gg(2)=yjk*fac/rjk
2163             gg(3)=zjk*fac/rjk
2164       do m=1,3
2165         gvdwx(m,j)=gvdwx(m,j)-gg(m)
2166         gvdwx(m,k)=gvdwx(m,k)+gg(m)
2167       enddo
2168       do l=1,3
2169         gvdwc(l,j)=gvdwc(l,j)-gg(l)
2170         gvdwc(l,k)=gvdwc(l,k)+gg(l)
2171       enddo
2172       return
2173       end