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