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