wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-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       include 'COMMON.MD'
102 #endif
103 #endif
104
105 c     External functions
106       double precision h_base
107       external h_base
108
109 c     Input arguments
110       integer resi,resj
111
112 c     Output arguments
113       double precision eij
114
115 c     Local variables
116       logical havebond
117 c      integer itypi,itypj,k,l
118       double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
119       double precision sig0ij,ljd,sig,fac,e1,e2
120       double precision dcosom1(3),dcosom2(3),ed
121       double precision pom1,pom2
122       double precision ljA,ljB,ljXs
123       double precision d_ljB(1:3)
124       double precision ssA,ssB,ssC,ssXs
125       double precision ssxm,ljxm,ssm,ljm
126       double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
127       double precision f1,f2,h1,h2,hd1,hd2
128       double precision omega,delta_inv,deltasq_inv,fac1,fac2
129 c-------FIRST METHOD
130       double precision xm,d_xm(1:3)
131 c-------END FIRST METHOD
132 c-------SECOND METHOD
133 c$$$      double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
134 c-------END SECOND METHOD
135
136 c-------TESTING CODE
137       logical checkstop,transgrad
138       common /sschecks/ checkstop,transgrad
139
140       integer icheck,nicheck,jcheck,njcheck
141       double precision echeck(-1:1),deps,ssx0,ljx0
142 c-------END TESTING CODE
143
144
145       i=resi
146       j=resj
147
148       itypi=itype(i)
149       dxi=dc_norm(1,nres+i)
150       dyi=dc_norm(2,nres+i)
151       dzi=dc_norm(3,nres+i)
152       dsci_inv=vbld_inv(i+nres)
153         xi=c(1,nres+i)
154         yi=c(2,nres+i)
155         zi=c(3,nres+i)
156           xi=mod(xi,boxxsize)
157           if (xi.lt.0) xi=xi+boxxsize
158           yi=mod(yi,boxysize)
159           if (yi.lt.0) yi=yi+boxysize
160           zi=mod(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=mod(xj,boxxsize)
195           if (xj.lt.0) xj=xj+boxxsize
196           yj=mod(yj,boxysize)
197           if (yj.lt.0) yj=yj+boxysize
198           zj=mod(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         if (.not.found.and.fg_rank.eq.0) 
746      &      write(iout,'(a15,f12.2,f8.1,2i5)')
747      &       "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         if (.not.found.and.fg_rank.eq.0) 
761      &      write(iout,'(a15,f12.2,f8.1,2i5)')
762      &       "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 c----------------------------------------------------------------------------
777
778 #ifdef WHAM
779       subroutine read_ssHist
780       implicit none
781
782 c     Includes
783       include 'DIMENSIONS'
784       include "DIMENSIONS.FREE"
785       include 'COMMON.FREE'
786
787 c     Local variables
788       integer i,j
789       character*80 controlcard
790
791       do i=1,dyn_nssHist
792         call card_concat(controlcard,.true.)
793         read(controlcard,*)
794      &       dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
795       enddo
796
797       return
798       end
799 #endif
800
801 c----------------------------------------------------------------------------
802
803
804 C-----------------------------------------------------------------------------
805 C-----------------------------------------------------------------------------
806 C-----------------------------------------------------------------------------
807 C-----------------------------------------------------------------------------
808 C-----------------------------------------------------------------------------
809 C-----------------------------------------------------------------------------
810 C-----------------------------------------------------------------------------
811
812 c$$$c-----------------------------------------------------------------------------
813 c$$$
814 c$$$      subroutine ss_relax(i_in,j_in)
815 c$$$      implicit none
816 c$$$
817 c$$$c     Includes
818 c$$$      include 'DIMENSIONS'
819 c$$$      include 'COMMON.VAR'
820 c$$$      include 'COMMON.CHAIN'
821 c$$$      include 'COMMON.IOUNITS'
822 c$$$      include 'COMMON.INTERACT'
823 c$$$
824 c$$$c     Input arguments
825 c$$$      integer i_in,j_in
826 c$$$
827 c$$$c     Local variables
828 c$$$      integer i,iretcode,nfun_sc
829 c$$$      logical scfail
830 c$$$      double precision var(maxvar),e_sc,etot
831 c$$$
832 c$$$
833 c$$$      mask_r=.true.
834 c$$$      do i=nnt,nct
835 c$$$        mask_side(i)=0
836 c$$$      enddo
837 c$$$      mask_side(i_in)=1
838 c$$$      mask_side(j_in)=1
839 c$$$
840 c$$$c     Minimize the two selected side-chains
841 c$$$      call overlap_sc(scfail)  ! Better not fail!
842 c$$$      call minimize_sc(e_sc,var,iretcode,nfun_sc)
843 c$$$
844 c$$$      mask_r=.false.
845 c$$$
846 c$$$      return
847 c$$$      end
848 c$$$
849 c$$$c-------------------------------------------------------------
850 c$$$
851 c$$$      subroutine minimize_sc(etot_sc,iretcode,nfun)
852 c$$$c     Minimize side-chains only, starting from geom but without modifying
853 c$$$c     bond lengths.
854 c$$$c     If mask_r is already set, only the selected side-chains are minimized,
855 c$$$c     otherwise all side-chains are minimized keeping the backbone frozen.
856 c$$$      implicit none
857 c$$$
858 c$$$c     Includes
859 c$$$      include 'DIMENSIONS'
860 c$$$      include 'COMMON.IOUNITS'
861 c$$$      include 'COMMON.VAR'
862 c$$$      include 'COMMON.CHAIN'
863 c$$$      include 'COMMON.GEO'
864 c$$$      include 'COMMON.MINIM'
865 c$$$      integer icall
866 c$$$      common /srutu/ icall
867 c$$$
868 c$$$c     Output arguments
869 c$$$      double precision etot_sc
870 c$$$      integer iretcode,nfun
871 c$$$
872 c$$$c     External functions/subroutines
873 c$$$      external func_sc,grad_sc,fdum
874 c$$$
875 c$$$c     Local variables
876 c$$$      integer liv,lv
877 c$$$      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
878 c$$$      integer iv(liv)
879 c$$$      double precision rdum(1)
880 c$$$      double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
881 c$$$      integer idum(1)
882 c$$$      integer i,nvar_restr
883 c$$$
884 c$$$
885 c$$$cmc      start_minim=.true.
886 c$$$      call deflt(2,iv,liv,lv,v)                                         
887 c$$$* 12 means fresh start, dont call deflt                                 
888 c$$$      iv(1)=12                                                          
889 c$$$* max num of fun calls                                                  
890 c$$$      if (maxfun.eq.0) maxfun=500
891 c$$$      iv(17)=maxfun
892 c$$$* max num of iterations                                                 
893 c$$$      if (maxmin.eq.0) maxmin=1000
894 c$$$      iv(18)=maxmin
895 c$$$* controls output                                                       
896 c$$$      iv(19)=1
897 c$$$* selects output unit                                                   
898 c$$$      iv(21)=0
899 c$$$c      iv(21)=iout               ! DEBUG
900 c$$$c      iv(21)=8                  ! DEBUG
901 c$$$* 1 means to print out result                                           
902 c$$$      iv(22)=0
903 c$$$c      iv(22)=1                  ! DEBUG
904 c$$$* 1 means to print out summary stats                                    
905 c$$$      iv(23)=0                                                          
906 c$$$c      iv(23)=1                  ! DEBUG
907 c$$$* 1 means to print initial x and d                                      
908 c$$$      iv(24)=0                                                          
909 c$$$c      iv(24)=1                  ! DEBUG
910 c$$$* min val for v(radfac) default is 0.1                                  
911 c$$$      v(24)=0.1D0                                                       
912 c$$$* max val for v(radfac) default is 4.0                                  
913 c$$$      v(25)=2.0D0                                                       
914 c$$$c     v(25)=4.0D0                                                       
915 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
916 c$$$* the sumsl default is 0.1                                              
917 c$$$      v(26)=0.1D0
918 c$$$* false conv if (act fnctn decrease) .lt. v(34)                         
919 c$$$* the sumsl default is 100*machep                                       
920 c$$$      v(34)=v(34)/100.0D0                                               
921 c$$$* absolute convergence                                                  
922 c$$$      if (tolf.eq.0.0D0) tolf=1.0D-4
923 c$$$      v(31)=tolf
924 c$$$* relative convergence                                                  
925 c$$$      if (rtolf.eq.0.0D0) rtolf=1.0D-1
926 c$$$      v(32)=rtolf
927 c$$$* controls initial step size                                            
928 c$$$       v(35)=1.0D-1                                                    
929 c$$$* large vals of d correspond to small components of step                
930 c$$$      do i=1,nphi
931 c$$$        d(i)=1.0D-1
932 c$$$      enddo
933 c$$$      do i=nphi+1,nvar
934 c$$$        d(i)=1.0D-1
935 c$$$      enddo
936 c$$$
937 c$$$      call geom_to_var(nvar,x)
938 c$$$      IF (mask_r) THEN
939 c$$$        do i=1,nres             ! Just in case...
940 c$$$          mask_phi(i)=0
941 c$$$          mask_theta(i)=0
942 c$$$        enddo
943 c$$$        call x2xx(x,xx,nvar_restr)
944 c$$$        call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
945 c$$$     &       iv,liv,lv,v,idum,rdum,fdum)      
946 c$$$        call xx2x(x,xx)
947 c$$$      ELSE
948 c$$$c     When minimizing ALL side-chains, etotal_sc is a little
949 c$$$c     faster if we don't set mask_r
950 c$$$        do i=1,nres
951 c$$$          mask_phi(i)=0
952 c$$$          mask_theta(i)=0
953 c$$$          mask_side(i)=1
954 c$$$        enddo
955 c$$$        call x2xx(x,xx,nvar_restr)
956 c$$$        call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
957 c$$$     &       iv,liv,lv,v,idum,rdum,fdum)      
958 c$$$        call xx2x(x,xx)
959 c$$$      ENDIF
960 c$$$      call var_to_geom(nvar,x)
961 c$$$      call chainbuild_sc
962 c$$$      etot_sc=v(10)                                                      
963 c$$$      iretcode=iv(1)
964 c$$$      nfun=iv(6)
965 c$$$      return  
966 c$$$      end  
967 c$$$
968 c$$$C--------------------------------------------------------------------------
969 c$$$
970 c$$$      subroutine chainbuild_sc
971 c$$$      implicit none
972 c$$$      include 'DIMENSIONS'
973 c$$$      include 'COMMON.VAR'
974 c$$$      include 'COMMON.INTERACT'
975 c$$$
976 c$$$c     Local variables
977 c$$$      integer i
978 c$$$
979 c$$$
980 c$$$      do i=nnt,nct
981 c$$$        if (.not.mask_r .or. mask_side(i).eq.1) then
982 c$$$          call locate_side_chain(i)
983 c$$$        endif
984 c$$$      enddo
985 c$$$
986 c$$$      return
987 c$$$      end
988 c$$$
989 c$$$C--------------------------------------------------------------------------
990 c$$$
991 c$$$      subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)  
992 c$$$      implicit none
993 c$$$
994 c$$$c     Includes
995 c$$$      include 'DIMENSIONS'
996 c$$$      include 'COMMON.DERIV'
997 c$$$      include 'COMMON.VAR'
998 c$$$      include 'COMMON.MINIM'
999 c$$$      include 'COMMON.IOUNITS'
1000 c$$$
1001 c$$$c     Input arguments
1002 c$$$      integer n
1003 c$$$      double precision x(maxvar)
1004 c$$$      double precision ufparm
1005 c$$$      external ufparm
1006 c$$$
1007 c$$$c     Input/Output arguments
1008 c$$$      integer nf
1009 c$$$      integer uiparm(1)
1010 c$$$      double precision urparm(1)
1011 c$$$
1012 c$$$c     Output arguments
1013 c$$$      double precision f
1014 c$$$
1015 c$$$c     Local variables
1016 c$$$      double precision energia(0:n_ene)
1017 c$$$#ifdef OSF
1018 c$$$c     Variables used to intercept NaNs
1019 c$$$      double precision x_sum
1020 c$$$      integer i_NAN
1021 c$$$#endif
1022 c$$$
1023 c$$$
1024 c$$$      nfl=nf
1025 c$$$      icg=mod(nf,2)+1
1026 c$$$
1027 c$$$#ifdef OSF
1028 c$$$c     Intercept NaNs in the coordinates, before calling etotal_sc
1029 c$$$      x_sum=0.D0
1030 c$$$      do i_NAN=1,n
1031 c$$$        x_sum=x_sum+x(i_NAN)
1032 c$$$      enddo
1033 c$$$c     Calculate the energy only if the coordinates are ok
1034 c$$$      if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
1035 c$$$        write(iout,*)"   *** func_restr_sc : Found NaN in coordinates"
1036 c$$$        f=1.0D+77
1037 c$$$        nf=0
1038 c$$$      else
1039 c$$$#endif
1040 c$$$
1041 c$$$      call var_to_geom_restr(n,x)
1042 c$$$      call zerograd
1043 c$$$      call chainbuild_sc
1044 c$$$      call etotal_sc(energia(0))
1045 c$$$      f=energia(0)
1046 c$$$      if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
1047 c$$$
1048 c$$$#ifdef OSF
1049 c$$$      endif
1050 c$$$#endif
1051 c$$$
1052 c$$$      return                                                            
1053 c$$$      end                                                               
1054 c$$$
1055 c$$$c-------------------------------------------------------
1056 c$$$
1057 c$$$      subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
1058 c$$$      implicit none
1059 c$$$
1060 c$$$c     Includes
1061 c$$$      include 'DIMENSIONS'
1062 c$$$      include 'COMMON.CHAIN'
1063 c$$$      include 'COMMON.DERIV'
1064 c$$$      include 'COMMON.VAR'
1065 c$$$      include 'COMMON.INTERACT'
1066 c$$$      include 'COMMON.MINIM'
1067 c$$$
1068 c$$$c     Input arguments
1069 c$$$      integer n
1070 c$$$      double precision x(maxvar)
1071 c$$$      double precision ufparm
1072 c$$$      external ufparm
1073 c$$$
1074 c$$$c     Input/Output arguments
1075 c$$$      integer nf
1076 c$$$      integer uiparm(1)
1077 c$$$      double precision urparm(1)
1078 c$$$
1079 c$$$c     Output arguments
1080 c$$$      double precision g(maxvar)
1081 c$$$
1082 c$$$c     Local variables
1083 c$$$      double precision f,gphii,gthetai,galphai,gomegai
1084 c$$$      integer ig,ind,i,j,k,igall,ij
1085 c$$$
1086 c$$$
1087 c$$$      icg=mod(nf,2)+1
1088 c$$$      if (nf-nfl+1) 20,30,40
1089 c$$$   20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
1090 c$$$c     write (iout,*) 'grad 20'
1091 c$$$      if (nf.eq.0) return
1092 c$$$      goto 40
1093 c$$$   30 call var_to_geom_restr(n,x)
1094 c$$$      call chainbuild_sc
1095 c$$$C
1096 c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
1097 c$$$C
1098 c$$$   40 call cartder
1099 c$$$C
1100 c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
1101 c$$$C
1102 c$$$
1103 c$$$      ig=0
1104 c$$$      ind=nres-2
1105 c$$$      do i=2,nres-2
1106 c$$$       IF (mask_phi(i+2).eq.1) THEN
1107 c$$$        gphii=0.0D0
1108 c$$$        do j=i+1,nres-1
1109 c$$$          ind=ind+1
1110 c$$$          do k=1,3
1111 c$$$            gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
1112 c$$$            gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
1113 c$$$          enddo
1114 c$$$        enddo
1115 c$$$        ig=ig+1
1116 c$$$        g(ig)=gphii
1117 c$$$       ELSE
1118 c$$$        ind=ind+nres-1-i
1119 c$$$       ENDIF
1120 c$$$      enddo                                        
1121 c$$$
1122 c$$$
1123 c$$$      ind=0
1124 c$$$      do i=1,nres-2
1125 c$$$       IF (mask_theta(i+2).eq.1) THEN
1126 c$$$        ig=ig+1
1127 c$$$    gthetai=0.0D0
1128 c$$$    do j=i+1,nres-1
1129 c$$$          ind=ind+1
1130 c$$$      do k=1,3
1131 c$$$            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
1132 c$$$            gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
1133 c$$$          enddo
1134 c$$$        enddo
1135 c$$$        g(ig)=gthetai
1136 c$$$       ELSE
1137 c$$$        ind=ind+nres-1-i
1138 c$$$       ENDIF
1139 c$$$      enddo
1140 c$$$
1141 c$$$      do i=2,nres-1
1142 c$$$    if (itype(i).ne.10) then
1143 c$$$         IF (mask_side(i).eq.1) THEN
1144 c$$$          ig=ig+1
1145 c$$$          galphai=0.0D0
1146 c$$$      do k=1,3
1147 c$$$        galphai=galphai+dxds(k,i)*gradx(k,i,icg)
1148 c$$$          enddo
1149 c$$$          g(ig)=galphai
1150 c$$$         ENDIF
1151 c$$$        endif
1152 c$$$      enddo
1153 c$$$
1154 c$$$      
1155 c$$$      do i=2,nres-1
1156 c$$$        if (itype(i).ne.10) then
1157 c$$$         IF (mask_side(i).eq.1) THEN
1158 c$$$          ig=ig+1
1159 c$$$      gomegai=0.0D0
1160 c$$$      do k=1,3
1161 c$$$        gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
1162 c$$$          enddo
1163 c$$$      g(ig)=gomegai
1164 c$$$         ENDIF
1165 c$$$        endif
1166 c$$$      enddo
1167 c$$$
1168 c$$$C
1169 c$$$C Add the components corresponding to local energy terms.
1170 c$$$C
1171 c$$$
1172 c$$$      ig=0
1173 c$$$      igall=0
1174 c$$$      do i=4,nres
1175 c$$$        igall=igall+1
1176 c$$$        if (mask_phi(i).eq.1) then
1177 c$$$          ig=ig+1
1178 c$$$          g(ig)=g(ig)+gloc(igall,icg)
1179 c$$$        endif
1180 c$$$      enddo
1181 c$$$
1182 c$$$      do i=3,nres
1183 c$$$        igall=igall+1
1184 c$$$        if (mask_theta(i).eq.1) then
1185 c$$$          ig=ig+1
1186 c$$$          g(ig)=g(ig)+gloc(igall,icg)
1187 c$$$        endif
1188 c$$$      enddo
1189 c$$$     
1190 c$$$      do ij=1,2
1191 c$$$      do i=2,nres-1
1192 c$$$        if (itype(i).ne.10) then
1193 c$$$          igall=igall+1
1194 c$$$          if (mask_side(i).eq.1) then
1195 c$$$            ig=ig+1
1196 c$$$            g(ig)=g(ig)+gloc(igall,icg)
1197 c$$$          endif
1198 c$$$        endif
1199 c$$$      enddo
1200 c$$$      enddo
1201 c$$$
1202 c$$$cd      do i=1,ig
1203 c$$$cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
1204 c$$$cd      enddo
1205 c$$$
1206 c$$$      return
1207 c$$$      end
1208 c$$$
1209 c$$$C-----------------------------------------------------------------------------
1210 c$$$
1211 c$$$      subroutine etotal_sc(energy_sc)
1212 c$$$      implicit none
1213 c$$$
1214 c$$$c     Includes
1215 c$$$      include 'DIMENSIONS'
1216 c$$$      include 'COMMON.VAR'
1217 c$$$      include 'COMMON.INTERACT'
1218 c$$$      include 'COMMON.DERIV'
1219 c$$$      include 'COMMON.FFIELD'
1220 c$$$
1221 c$$$c     Output arguments
1222 c$$$      double precision energy_sc(0:n_ene)
1223 c$$$
1224 c$$$c     Local variables
1225 c$$$      double precision evdw,escloc
1226 c$$$      integer i,j
1227 c$$$
1228 c$$$
1229 c$$$      do i=1,n_ene
1230 c$$$        energy_sc(i)=0.0D0
1231 c$$$      enddo
1232 c$$$
1233 c$$$      if (mask_r) then
1234 c$$$        call egb_sc(evdw)
1235 c$$$        call esc_sc(escloc)
1236 c$$$      else
1237 c$$$        call egb(evdw)
1238 c$$$        call esc(escloc)
1239 c$$$      endif
1240 c$$$
1241 c$$$      if (evdw.eq.1.0D20) then
1242 c$$$        energy_sc(0)=evdw
1243 c$$$      else
1244 c$$$        energy_sc(0)=wsc*evdw+wscloc*escloc
1245 c$$$      endif
1246 c$$$      energy_sc(1)=evdw
1247 c$$$      energy_sc(12)=escloc
1248 c$$$
1249 c$$$C
1250 c$$$C Sum up the components of the Cartesian gradient.
1251 c$$$C
1252 c$$$      do i=1,nct
1253 c$$$        do j=1,3
1254 c$$$          gradx(j,i,icg)=wsc*gvdwx(j,i)
1255 c$$$        enddo
1256 c$$$      enddo
1257 c$$$
1258 c$$$      return
1259 c$$$      end
1260 c$$$
1261 c$$$C-----------------------------------------------------------------------------
1262 c$$$
1263 c$$$      subroutine egb_sc(evdw)
1264 c$$$C
1265 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1266 c$$$C assuming the Gay-Berne potential of interaction.
1267 c$$$C
1268 c$$$      implicit real*8 (a-h,o-z)
1269 c$$$      include 'DIMENSIONS'
1270 c$$$      include 'COMMON.GEO'
1271 c$$$      include 'COMMON.VAR'
1272 c$$$      include 'COMMON.LOCAL'
1273 c$$$      include 'COMMON.CHAIN'
1274 c$$$      include 'COMMON.DERIV'
1275 c$$$      include 'COMMON.NAMES'
1276 c$$$      include 'COMMON.INTERACT'
1277 c$$$      include 'COMMON.IOUNITS'
1278 c$$$      include 'COMMON.CALC'
1279 c$$$      include 'COMMON.CONTROL'
1280 c$$$      logical lprn
1281 c$$$      evdw=0.0D0
1282 c$$$      energy_dec=.false.
1283 c$$$c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1284 c$$$      evdw=0.0D0
1285 c$$$      lprn=.false.
1286 c$$$c     if (icall.eq.0) lprn=.false.
1287 c$$$      ind=0
1288 c$$$      do i=iatsc_s,iatsc_e
1289 c$$$        itypi=itype(i)
1290 c$$$        itypi1=itype(i+1)
1291 c$$$        xi=c(1,nres+i)
1292 c$$$        yi=c(2,nres+i)
1293 c$$$        zi=c(3,nres+i)
1294 c$$$        dxi=dc_norm(1,nres+i)
1295 c$$$        dyi=dc_norm(2,nres+i)
1296 c$$$        dzi=dc_norm(3,nres+i)
1297 c$$$c        dsci_inv=dsc_inv(itypi)
1298 c$$$        dsci_inv=vbld_inv(i+nres)
1299 c$$$c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1300 c$$$c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1301 c$$$C
1302 c$$$C Calculate SC interaction energy.
1303 c$$$C
1304 c$$$        do iint=1,nint_gr(i)
1305 c$$$          do j=istart(i,iint),iend(i,iint)
1306 c$$$          IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
1307 c$$$            ind=ind+1
1308 c$$$            itypj=itype(j)
1309 c$$$c            dscj_inv=dsc_inv(itypj)
1310 c$$$            dscj_inv=vbld_inv(j+nres)
1311 c$$$c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1312 c$$$c     &       1.0d0/vbld(j+nres)
1313 c$$$c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1314 c$$$            sig0ij=sigma(itypi,itypj)
1315 c$$$            chi1=chi(itypi,itypj)
1316 c$$$            chi2=chi(itypj,itypi)
1317 c$$$            chi12=chi1*chi2
1318 c$$$            chip1=chip(itypi)
1319 c$$$            chip2=chip(itypj)
1320 c$$$            chip12=chip1*chip2
1321 c$$$            alf1=alp(itypi)
1322 c$$$            alf2=alp(itypj)
1323 c$$$            alf12=0.5D0*(alf1+alf2)
1324 c$$$C For diagnostics only!!!
1325 c$$$c           chi1=0.0D0
1326 c$$$c           chi2=0.0D0
1327 c$$$c           chi12=0.0D0
1328 c$$$c           chip1=0.0D0
1329 c$$$c           chip2=0.0D0
1330 c$$$c           chip12=0.0D0
1331 c$$$c           alf1=0.0D0
1332 c$$$c           alf2=0.0D0
1333 c$$$c           alf12=0.0D0
1334 c$$$            xj=c(1,nres+j)-xi
1335 c$$$            yj=c(2,nres+j)-yi
1336 c$$$            zj=c(3,nres+j)-zi
1337 c$$$            dxj=dc_norm(1,nres+j)
1338 c$$$            dyj=dc_norm(2,nres+j)
1339 c$$$            dzj=dc_norm(3,nres+j)
1340 c$$$c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1341 c$$$c            write (iout,*) "j",j," dc_norm",
1342 c$$$c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1343 c$$$            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1344 c$$$            rij=dsqrt(rrij)
1345 c$$$C Calculate angle-dependent terms of energy and contributions to their
1346 c$$$C derivatives.
1347 c$$$            call sc_angular
1348 c$$$            sigsq=1.0D0/sigsq
1349 c$$$            sig=sig0ij*dsqrt(sigsq)
1350 c$$$            rij_shift=1.0D0/rij-sig+sig0ij
1351 c$$$c for diagnostics; uncomment
1352 c$$$c            rij_shift=1.2*sig0ij
1353 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1354 c$$$            if (rij_shift.le.0.0D0) then
1355 c$$$              evdw=1.0D20
1356 c$$$cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1357 c$$$cd     &        restyp(itypi),i,restyp(itypj),j,
1358 c$$$cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
1359 c$$$              return
1360 c$$$            endif
1361 c$$$            sigder=-sig*sigsq
1362 c$$$c---------------------------------------------------------------
1363 c$$$            rij_shift=1.0D0/rij_shift 
1364 c$$$            fac=rij_shift**expon
1365 c$$$            e1=fac*fac*aa(itypi,itypj)
1366 c$$$            e2=fac*bb(itypi,itypj)
1367 c$$$            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1368 c$$$            eps2der=evdwij*eps3rt
1369 c$$$            eps3der=evdwij*eps2rt
1370 c$$$c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1371 c$$$c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1372 c$$$            evdwij=evdwij*eps2rt*eps3rt
1373 c$$$            evdw=evdw+evdwij
1374 c$$$            if (lprn) then
1375 c$$$            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1376 c$$$            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1377 c$$$            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1378 c$$$     &        restyp(itypi),i,restyp(itypj),j,
1379 c$$$     &        epsi,sigm,chi1,chi2,chip1,chip2,
1380 c$$$     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1381 c$$$     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1382 c$$$     &        evdwij
1383 c$$$            endif
1384 c$$$
1385 c$$$            if (energy_dec) write (iout,'(a6,2i,0pf7.3)') 
1386 c$$$     &                        'evdw',i,j,evdwij
1387 c$$$
1388 c$$$C Calculate gradient components.
1389 c$$$            e1=e1*eps1*eps2rt**2*eps3rt**2
1390 c$$$            fac=-expon*(e1+evdwij)*rij_shift
1391 c$$$            sigder=fac*sigder
1392 c$$$            fac=rij*fac
1393 c$$$c            fac=0.0d0
1394 c$$$C Calculate the radial part of the gradient
1395 c$$$            gg(1)=xj*fac
1396 c$$$            gg(2)=yj*fac
1397 c$$$            gg(3)=zj*fac
1398 c$$$C Calculate angular part of the gradient.
1399 c$$$            call sc_grad
1400 c$$$          ENDIF
1401 c$$$          enddo      ! j
1402 c$$$        enddo        ! iint
1403 c$$$      enddo          ! i
1404 c$$$      energy_dec=.false.
1405 c$$$      return
1406 c$$$      end
1407 c$$$
1408 c$$$c-----------------------------------------------------------------------------
1409 c$$$
1410 c$$$      subroutine esc_sc(escloc)
1411 c$$$C Calculate the local energy of a side chain and its derivatives in the
1412 c$$$C corresponding virtual-bond valence angles THETA and the spherical angles 
1413 c$$$C ALPHA and OMEGA.
1414 c$$$      implicit real*8 (a-h,o-z)
1415 c$$$      include 'DIMENSIONS'
1416 c$$$      include 'COMMON.GEO'
1417 c$$$      include 'COMMON.LOCAL'
1418 c$$$      include 'COMMON.VAR'
1419 c$$$      include 'COMMON.INTERACT'
1420 c$$$      include 'COMMON.DERIV'
1421 c$$$      include 'COMMON.CHAIN'
1422 c$$$      include 'COMMON.IOUNITS'
1423 c$$$      include 'COMMON.NAMES'
1424 c$$$      include 'COMMON.FFIELD'
1425 c$$$      include 'COMMON.CONTROL'
1426 c$$$      double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
1427 c$$$     &     ddersc0(3),ddummy(3),xtemp(3),temp(3)
1428 c$$$      common /sccalc/ time11,time12,time112,theti,it,nlobit
1429 c$$$      delta=0.02d0*pi
1430 c$$$      escloc=0.0D0
1431 c$$$c     write (iout,'(a)') 'ESC'
1432 c$$$      do i=loc_start,loc_end
1433 c$$$      IF (mask_side(i).eq.1) THEN
1434 c$$$        it=itype(i)
1435 c$$$        if (it.eq.10) goto 1
1436 c$$$        nlobit=nlob(it)
1437 c$$$c       print *,'i=',i,' it=',it,' nlobit=',nlobit
1438 c$$$c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
1439 c$$$        theti=theta(i+1)-pipol
1440 c$$$        x(1)=dtan(theti)
1441 c$$$        x(2)=alph(i)
1442 c$$$        x(3)=omeg(i)
1443 c$$$
1444 c$$$        if (x(2).gt.pi-delta) then
1445 c$$$          xtemp(1)=x(1)
1446 c$$$          xtemp(2)=pi-delta
1447 c$$$          xtemp(3)=x(3)
1448 c$$$          call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1449 c$$$          xtemp(2)=pi
1450 c$$$          call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1451 c$$$          call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
1452 c$$$     &        escloci,dersc(2))
1453 c$$$          call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1454 c$$$     &        ddersc0(1),dersc(1))
1455 c$$$          call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
1456 c$$$     &        ddersc0(3),dersc(3))
1457 c$$$          xtemp(2)=pi-delta
1458 c$$$          call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1459 c$$$          xtemp(2)=pi
1460 c$$$          call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1461 c$$$          call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
1462 c$$$     &            dersc0(2),esclocbi,dersc02)
1463 c$$$          call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
1464 c$$$     &            dersc12,dersc01)
1465 c$$$          call splinthet(x(2),0.5d0*delta,ss,ssd)
1466 c$$$          dersc0(1)=dersc01
1467 c$$$          dersc0(2)=dersc02
1468 c$$$          dersc0(3)=0.0d0
1469 c$$$          do k=1,3
1470 c$$$            dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1471 c$$$          enddo
1472 c$$$          dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1473 c$$$c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1474 c$$$c    &             esclocbi,ss,ssd
1475 c$$$          escloci=ss*escloci+(1.0d0-ss)*esclocbi
1476 c$$$c         escloci=esclocbi
1477 c$$$c         write (iout,*) escloci
1478 c$$$        else if (x(2).lt.delta) then
1479 c$$$          xtemp(1)=x(1)
1480 c$$$          xtemp(2)=delta
1481 c$$$          xtemp(3)=x(3)
1482 c$$$          call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
1483 c$$$          xtemp(2)=0.0d0
1484 c$$$          call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
1485 c$$$          call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
1486 c$$$     &        escloci,dersc(2))
1487 c$$$          call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1488 c$$$     &        ddersc0(1),dersc(1))
1489 c$$$          call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
1490 c$$$     &        ddersc0(3),dersc(3))
1491 c$$$          xtemp(2)=delta
1492 c$$$          call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
1493 c$$$          xtemp(2)=0.0d0
1494 c$$$          call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
1495 c$$$          call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
1496 c$$$     &            dersc0(2),esclocbi,dersc02)
1497 c$$$          call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
1498 c$$$     &            dersc12,dersc01)
1499 c$$$          dersc0(1)=dersc01
1500 c$$$          dersc0(2)=dersc02
1501 c$$$          dersc0(3)=0.0d0
1502 c$$$          call splinthet(x(2),0.5d0*delta,ss,ssd)
1503 c$$$          do k=1,3
1504 c$$$            dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
1505 c$$$          enddo
1506 c$$$          dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
1507 c$$$c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
1508 c$$$c    &             esclocbi,ss,ssd
1509 c$$$          escloci=ss*escloci+(1.0d0-ss)*esclocbi
1510 c$$$c         write (iout,*) escloci
1511 c$$$        else
1512 c$$$          call enesc(x,escloci,dersc,ddummy,.false.)
1513 c$$$        endif
1514 c$$$
1515 c$$$        escloc=escloc+escloci
1516 c$$$        if (energy_dec) write (iout,'(a6,i,0pf7.3)')
1517 c$$$     &     'escloc',i,escloci
1518 c$$$c       write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
1519 c$$$
1520 c$$$        gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
1521 c$$$     &   wscloc*dersc(1)
1522 c$$$        gloc(ialph(i,1),icg)=wscloc*dersc(2)
1523 c$$$        gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
1524 c$$$    1   continue
1525 c$$$      ENDIF
1526 c$$$      enddo
1527 c$$$      return
1528 c$$$      end
1529 c$$$
1530 c$$$C-----------------------------------------------------------------------------
1531 c$$$
1532 c$$$      subroutine egb_ij(i_sc,j_sc,evdw)
1533 c$$$C
1534 c$$$C This subroutine calculates the interaction energy of nonbonded side chains
1535 c$$$C assuming the Gay-Berne potential of interaction.
1536 c$$$C
1537 c$$$      implicit real*8 (a-h,o-z)
1538 c$$$      include 'DIMENSIONS'
1539 c$$$      include 'COMMON.GEO'
1540 c$$$      include 'COMMON.VAR'
1541 c$$$      include 'COMMON.LOCAL'
1542 c$$$      include 'COMMON.CHAIN'
1543 c$$$      include 'COMMON.DERIV'
1544 c$$$      include 'COMMON.NAMES'
1545 c$$$      include 'COMMON.INTERACT'
1546 c$$$      include 'COMMON.IOUNITS'
1547 c$$$      include 'COMMON.CALC'
1548 c$$$      include 'COMMON.CONTROL'
1549 c$$$      logical lprn
1550 c$$$      evdw=0.0D0
1551 c$$$      energy_dec=.false.
1552 c$$$c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1553 c$$$      evdw=0.0D0
1554 c$$$      lprn=.false.
1555 c$$$      ind=0
1556 c$$$c$$$      do i=iatsc_s,iatsc_e
1557 c$$$      i=i_sc
1558 c$$$        itypi=itype(i)
1559 c$$$        itypi1=itype(i+1)
1560 c$$$        xi=c(1,nres+i)
1561 c$$$        yi=c(2,nres+i)
1562 c$$$        zi=c(3,nres+i)
1563 c$$$        dxi=dc_norm(1,nres+i)
1564 c$$$        dyi=dc_norm(2,nres+i)
1565 c$$$        dzi=dc_norm(3,nres+i)
1566 c$$$c        dsci_inv=dsc_inv(itypi)
1567 c$$$        dsci_inv=vbld_inv(i+nres)
1568 c$$$c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1569 c$$$c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1570 c$$$C
1571 c$$$C Calculate SC interaction energy.
1572 c$$$C
1573 c$$$c$$$        do iint=1,nint_gr(i)
1574 c$$$c$$$          do j=istart(i,iint),iend(i,iint)
1575 c$$$        j=j_sc
1576 c$$$            ind=ind+1
1577 c$$$            itypj=itype(j)
1578 c$$$c            dscj_inv=dsc_inv(itypj)
1579 c$$$            dscj_inv=vbld_inv(j+nres)
1580 c$$$c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1581 c$$$c     &       1.0d0/vbld(j+nres)
1582 c$$$c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1583 c$$$            sig0ij=sigma(itypi,itypj)
1584 c$$$            chi1=chi(itypi,itypj)
1585 c$$$            chi2=chi(itypj,itypi)
1586 c$$$            chi12=chi1*chi2
1587 c$$$            chip1=chip(itypi)
1588 c$$$            chip2=chip(itypj)
1589 c$$$            chip12=chip1*chip2
1590 c$$$            alf1=alp(itypi)
1591 c$$$            alf2=alp(itypj)
1592 c$$$            alf12=0.5D0*(alf1+alf2)
1593 c$$$C For diagnostics only!!!
1594 c$$$c           chi1=0.0D0
1595 c$$$c           chi2=0.0D0
1596 c$$$c           chi12=0.0D0
1597 c$$$c           chip1=0.0D0
1598 c$$$c           chip2=0.0D0
1599 c$$$c           chip12=0.0D0
1600 c$$$c           alf1=0.0D0
1601 c$$$c           alf2=0.0D0
1602 c$$$c           alf12=0.0D0
1603 c$$$            xj=c(1,nres+j)-xi
1604 c$$$            yj=c(2,nres+j)-yi
1605 c$$$            zj=c(3,nres+j)-zi
1606 c$$$            dxj=dc_norm(1,nres+j)
1607 c$$$            dyj=dc_norm(2,nres+j)
1608 c$$$            dzj=dc_norm(3,nres+j)
1609 c$$$c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
1610 c$$$c            write (iout,*) "j",j," dc_norm",
1611 c$$$c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
1612 c$$$            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1613 c$$$            rij=dsqrt(rrij)
1614 c$$$C Calculate angle-dependent terms of energy and contributions to their
1615 c$$$C derivatives.
1616 c$$$            call sc_angular
1617 c$$$            sigsq=1.0D0/sigsq
1618 c$$$            sig=sig0ij*dsqrt(sigsq)
1619 c$$$            rij_shift=1.0D0/rij-sig+sig0ij
1620 c$$$c for diagnostics; uncomment
1621 c$$$c            rij_shift=1.2*sig0ij
1622 c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
1623 c$$$            if (rij_shift.le.0.0D0) then
1624 c$$$              evdw=1.0D20
1625 c$$$cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1626 c$$$cd     &        restyp(itypi),i,restyp(itypj),j,
1627 c$$$cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
1628 c$$$              return
1629 c$$$            endif
1630 c$$$            sigder=-sig*sigsq
1631 c$$$c---------------------------------------------------------------
1632 c$$$            rij_shift=1.0D0/rij_shift 
1633 c$$$            fac=rij_shift**expon
1634 c$$$            e1=fac*fac*aa(itypi,itypj)
1635 c$$$            e2=fac*bb(itypi,itypj)
1636 c$$$            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1637 c$$$            eps2der=evdwij*eps3rt
1638 c$$$            eps3der=evdwij*eps2rt
1639 c$$$c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1640 c$$$c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1641 c$$$            evdwij=evdwij*eps2rt*eps3rt
1642 c$$$            evdw=evdw+evdwij
1643 c$$$            if (lprn) then
1644 c$$$            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1645 c$$$            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1646 c$$$            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1647 c$$$     &        restyp(itypi),i,restyp(itypj),j,
1648 c$$$     &        epsi,sigm,chi1,chi2,chip1,chip2,
1649 c$$$     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1650 c$$$     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1651 c$$$     &        evdwij
1652 c$$$            endif
1653 c$$$
1654 c$$$            if (energy_dec) write (iout,'(a6,2i,0pf7.3)') 
1655 c$$$     &                        'evdw',i,j,evdwij
1656 c$$$
1657 c$$$C Calculate gradient components.
1658 c$$$            e1=e1*eps1*eps2rt**2*eps3rt**2
1659 c$$$            fac=-expon*(e1+evdwij)*rij_shift
1660 c$$$            sigder=fac*sigder
1661 c$$$            fac=rij*fac
1662 c$$$c            fac=0.0d0
1663 c$$$C Calculate the radial part of the gradient
1664 c$$$            gg(1)=xj*fac
1665 c$$$            gg(2)=yj*fac
1666 c$$$            gg(3)=zj*fac
1667 c$$$C Calculate angular part of the gradient.
1668 c$$$            call sc_grad
1669 c$$$c$$$          enddo      ! j
1670 c$$$c$$$        enddo        ! iint
1671 c$$$c$$$      enddo          ! i
1672 c$$$      energy_dec=.false.
1673 c$$$      return
1674 c$$$      end
1675 c$$$
1676 c$$$C-----------------------------------------------------------------------------
1677 c$$$
1678 c$$$      subroutine perturb_side_chain(i,angle)
1679 c$$$      implicit none
1680 c$$$
1681 c$$$c     Includes
1682 c$$$      include 'DIMENSIONS'
1683 c$$$      include 'COMMON.CHAIN'
1684 c$$$      include 'COMMON.GEO'
1685 c$$$      include 'COMMON.VAR'
1686 c$$$      include 'COMMON.LOCAL'
1687 c$$$      include 'COMMON.IOUNITS'
1688 c$$$
1689 c$$$c     External functions
1690 c$$$      external ran_number
1691 c$$$      double precision ran_number
1692 c$$$
1693 c$$$c     Input arguments
1694 c$$$      integer i
1695 c$$$      double precision angle    ! In degrees
1696 c$$$
1697 c$$$c     Local variables
1698 c$$$      integer i_sc
1699 c$$$      double precision rad_ang,rand_v(3),length,cost,sint
1700 c$$$
1701 c$$$
1702 c$$$      i_sc=i+nres
1703 c$$$      rad_ang=angle*deg2rad
1704 c$$$
1705 c$$$      length=0.0
1706 c$$$      do while (length.lt.0.01)
1707 c$$$        rand_v(1)=ran_number(0.01D0,1.0D0)
1708 c$$$        rand_v(2)=ran_number(0.01D0,1.0D0)
1709 c$$$        rand_v(3)=ran_number(0.01D0,1.0D0)
1710 c$$$        length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
1711 c$$$     +       rand_v(3)*rand_v(3)
1712 c$$$        length=sqrt(length)
1713 c$$$        rand_v(1)=rand_v(1)/length
1714 c$$$        rand_v(2)=rand_v(2)/length
1715 c$$$        rand_v(3)=rand_v(3)/length
1716 c$$$        cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
1717 c$$$     +       rand_v(3)*dc_norm(3,i_sc)
1718 c$$$        length=1.0D0-cost*cost
1719 c$$$        if (length.lt.0.0D0) length=0.0D0
1720 c$$$        length=sqrt(length)
1721 c$$$        rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
1722 c$$$        rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
1723 c$$$        rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
1724 c$$$      enddo
1725 c$$$      rand_v(1)=rand_v(1)/length
1726 c$$$      rand_v(2)=rand_v(2)/length
1727 c$$$      rand_v(3)=rand_v(3)/length
1728 c$$$
1729 c$$$      cost=dcos(rad_ang)
1730 c$$$      sint=dsin(rad_ang)
1731 c$$$      dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
1732 c$$$      dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
1733 c$$$      dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
1734 c$$$      dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
1735 c$$$      dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
1736 c$$$      dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
1737 c$$$      c(1,i_sc)=c(1,i)+dc(1,i_sc)
1738 c$$$      c(2,i_sc)=c(2,i)+dc(2,i_sc)
1739 c$$$      c(3,i_sc)=c(3,i)+dc(3,i_sc)
1740 c$$$
1741 c$$$      call chainbuild_cart
1742 c$$$
1743 c$$$      return
1744 c$$$      end
1745 c$$$
1746 c$$$c----------------------------------------------------------------------------
1747 c$$$
1748 c$$$      subroutine ss_relax3(i_in,j_in)
1749 c$$$      implicit none
1750 c$$$
1751 c$$$c     Includes
1752 c$$$      include 'DIMENSIONS'
1753 c$$$      include 'COMMON.VAR'
1754 c$$$      include 'COMMON.CHAIN'
1755 c$$$      include 'COMMON.IOUNITS'
1756 c$$$      include 'COMMON.INTERACT'
1757 c$$$
1758 c$$$c     External functions
1759 c$$$      external ran_number
1760 c$$$      double precision ran_number
1761 c$$$
1762 c$$$c     Input arguments
1763 c$$$      integer i_in,j_in
1764 c$$$
1765 c$$$c     Local variables
1766 c$$$      double precision energy_sc(0:n_ene),etot
1767 c$$$      double precision org_dc(3),org_dc_norm(3),org_c(3)
1768 c$$$      double precision ang_pert,rand_fact,exp_fact,beta
1769 c$$$      integer n,i_pert,i
1770 c$$$      logical notdone
1771 c$$$
1772 c$$$
1773 c$$$      beta=1.0D0
1774 c$$$
1775 c$$$      mask_r=.true.
1776 c$$$      do i=nnt,nct
1777 c$$$        mask_side(i)=0
1778 c$$$      enddo
1779 c$$$      mask_side(i_in)=1
1780 c$$$      mask_side(j_in)=1
1781 c$$$
1782 c$$$      call etotal_sc(energy_sc)
1783 c$$$      etot=energy_sc(0)
1784 c$$$c      write(iout,'(a,3d15.5)')"     SS_MC_START ",energy_sc(0),
1785 c$$$c     +     energy_sc(1),energy_sc(12)
1786 c$$$
1787 c$$$      notdone=.true.
1788 c$$$      n=0
1789 c$$$      do while (notdone)
1790 c$$$        if (mod(n,2).eq.0) then
1791 c$$$          i_pert=i_in
1792 c$$$        else
1793 c$$$          i_pert=j_in
1794 c$$$        endif
1795 c$$$        n=n+1
1796 c$$$
1797 c$$$        do i=1,3
1798 c$$$          org_dc(i)=dc(i,i_pert+nres)
1799 c$$$          org_dc_norm(i)=dc_norm(i,i_pert+nres)
1800 c$$$          org_c(i)=c(i,i_pert+nres)
1801 c$$$        enddo
1802 c$$$        ang_pert=ran_number(0.0D0,3.0D0)
1803 c$$$        call perturb_side_chain(i_pert,ang_pert)
1804 c$$$        call etotal_sc(energy_sc)
1805 c$$$        exp_fact=exp(beta*(etot-energy_sc(0)))
1806 c$$$        rand_fact=ran_number(0.0D0,1.0D0)
1807 c$$$        if (rand_fact.lt.exp_fact) then
1808 c$$$c          write(iout,'(a,3d15.5)')"     SS_MC_ACCEPT ",energy_sc(0),
1809 c$$$c     +     energy_sc(1),energy_sc(12)
1810 c$$$          etot=energy_sc(0)
1811 c$$$        else
1812 c$$$c          write(iout,'(a,3d15.5)')"     SS_MC_REJECT ",energy_sc(0),
1813 c$$$c     +     energy_sc(1),energy_sc(12)
1814 c$$$          do i=1,3
1815 c$$$            dc(i,i_pert+nres)=org_dc(i)
1816 c$$$            dc_norm(i,i_pert+nres)=org_dc_norm(i)
1817 c$$$            c(i,i_pert+nres)=org_c(i)
1818 c$$$          enddo
1819 c$$$        endif
1820 c$$$
1821 c$$$        if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
1822 c$$$      enddo
1823 c$$$
1824 c$$$      mask_r=.false.
1825 c$$$
1826 c$$$      return
1827 c$$$      end
1828 c$$$
1829 c$$$c----------------------------------------------------------------------------
1830 c$$$
1831 c$$$      subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
1832 c$$$      implicit none
1833 c$$$      include 'DIMENSIONS'
1834 c$$$      integer liv,lv
1835 c$$$      parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2)) 
1836 c$$$*********************************************************************
1837 c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
1838 c$$$* the calling subprogram.                                           *     
1839 c$$$* when d(i)=1.0, then v(35) is the length of the initial step,      *     
1840 c$$$* calculated in the usual pythagorean way.                          *     
1841 c$$$* absolute convergence occurs when the function is within v(31) of  *     
1842 c$$$* zero. unless you know the minimum value in advance, abs convg     *     
1843 c$$$* is probably not useful.                                           *     
1844 c$$$* relative convergence is when the model predicts that the function *   
1845 c$$$* will decrease by less than v(32)*abs(fun).                        *   
1846 c$$$*********************************************************************
1847 c$$$      include 'COMMON.IOUNITS'
1848 c$$$      include 'COMMON.VAR'
1849 c$$$      include 'COMMON.GEO'
1850 c$$$      include 'COMMON.MINIM'
1851 c$$$      include 'COMMON.CHAIN'
1852 c$$$
1853 c$$$      double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1854 c$$$      common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1855 c$$$     +     orig_ss_dist(maxres2,maxres2)
1856 c$$$
1857 c$$$      double precision etot
1858 c$$$      integer iretcode,nfun,i_in,j_in
1859 c$$$
1860 c$$$      external dist
1861 c$$$      double precision dist
1862 c$$$      external ss_func,fdum
1863 c$$$      double precision ss_func,fdum
1864 c$$$
1865 c$$$      integer iv(liv),uiparm(2)
1866 c$$$      double precision v(lv),x(maxres6),d(maxres6),rdum
1867 c$$$      integer i,j,k
1868 c$$$
1869 c$$$
1870 c$$$      call deflt(2,iv,liv,lv,v)                                         
1871 c$$$* 12 means fresh start, dont call deflt                                 
1872 c$$$      iv(1)=12                                                          
1873 c$$$* max num of fun calls                                                  
1874 c$$$      if (maxfun.eq.0) maxfun=500
1875 c$$$      iv(17)=maxfun
1876 c$$$* max num of iterations                                                 
1877 c$$$      if (maxmin.eq.0) maxmin=1000
1878 c$$$      iv(18)=maxmin
1879 c$$$* controls output                                                       
1880 c$$$      iv(19)=2                                                          
1881 c$$$* selects output unit                                                   
1882 c$$$c      iv(21)=iout                                                       
1883 c$$$      iv(21)=0
1884 c$$$* 1 means to print out result                                           
1885 c$$$      iv(22)=0                                                          
1886 c$$$* 1 means to print out summary stats                                    
1887 c$$$      iv(23)=0                                                          
1888 c$$$* 1 means to print initial x and d                                      
1889 c$$$      iv(24)=0                                                          
1890 c$$$* min val for v(radfac) default is 0.1                                  
1891 c$$$      v(24)=0.1D0                                                       
1892 c$$$* max val for v(radfac) default is 4.0                                  
1893 c$$$      v(25)=2.0D0                                                       
1894 c$$$c     v(25)=4.0D0                                                       
1895 c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
1896 c$$$* the sumsl default is 0.1                                              
1897 c$$$      v(26)=0.1D0
1898 c$$$* false conv if (act fnctn decrease) .lt. v(34)                         
1899 c$$$* the sumsl default is 100*machep                                       
1900 c$$$      v(34)=v(34)/100.0D0                                               
1901 c$$$* absolute convergence                                                  
1902 c$$$      if (tolf.eq.0.0D0) tolf=1.0D-4
1903 c$$$      v(31)=tolf
1904 c$$$      v(31)=1.0D-1
1905 c$$$* relative convergence                                                  
1906 c$$$      if (rtolf.eq.0.0D0) rtolf=1.0D-4
1907 c$$$      v(32)=rtolf
1908 c$$$      v(32)=1.0D-1
1909 c$$$* controls initial step size                                            
1910 c$$$      v(35)=1.0D-1
1911 c$$$* large vals of d correspond to small components of step                
1912 c$$$      do i=1,6*nres
1913 c$$$        d(i)=1.0D0
1914 c$$$      enddo
1915 c$$$
1916 c$$$      do i=0,2*nres
1917 c$$$        do j=1,3
1918 c$$$          orig_ss_dc(j,i)=dc(j,i)
1919 c$$$        enddo
1920 c$$$      enddo
1921 c$$$      call geom_to_var(nvar,orig_ss_var)
1922 c$$$
1923 c$$$      do i=1,nres
1924 c$$$        do j=i,nres
1925 c$$$          orig_ss_dist(j,i)=dist(j,i)
1926 c$$$          orig_ss_dist(j+nres,i)=dist(j+nres,i)
1927 c$$$          orig_ss_dist(j,i+nres)=dist(j,i+nres)
1928 c$$$          orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
1929 c$$$        enddo
1930 c$$$      enddo
1931 c$$$
1932 c$$$      k=0
1933 c$$$      do i=1,nres-1
1934 c$$$        do j=1,3
1935 c$$$          k=k+1
1936 c$$$          x(k)=dc(j,i)
1937 c$$$        enddo
1938 c$$$      enddo
1939 c$$$      do i=2,nres-1
1940 c$$$        if (ialph(i,1).gt.0) then
1941 c$$$        do j=1,3
1942 c$$$          k=k+1
1943 c$$$          x(k)=dc(j,i+nres)
1944 c$$$        enddo
1945 c$$$        endif
1946 c$$$      enddo
1947 c$$$
1948 c$$$      uiparm(1)=i_in
1949 c$$$      uiparm(2)=j_in
1950 c$$$      call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
1951 c$$$      etot=v(10)
1952 c$$$      iretcode=iv(1)
1953 c$$$      nfun=iv(6)+iv(30)
1954 c$$$
1955 c$$$      k=0
1956 c$$$      do i=1,nres-1
1957 c$$$        do j=1,3
1958 c$$$          k=k+1
1959 c$$$          dc(j,i)=x(k)
1960 c$$$        enddo
1961 c$$$      enddo
1962 c$$$      do i=2,nres-1
1963 c$$$        if (ialph(i,1).gt.0) then
1964 c$$$        do j=1,3
1965 c$$$          k=k+1
1966 c$$$          dc(j,i+nres)=x(k)
1967 c$$$        enddo
1968 c$$$        endif
1969 c$$$      enddo
1970 c$$$      call chainbuild_cart
1971 c$$$
1972 c$$$      return  
1973 c$$$      end  
1974 c$$$
1975 c$$$C-----------------------------------------------------------------------------
1976 c$$$
1977 c$$$      subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)  
1978 c$$$      implicit none
1979 c$$$      include 'DIMENSIONS'
1980 c$$$      include 'COMMON.DERIV'
1981 c$$$      include 'COMMON.IOUNITS'
1982 c$$$      include 'COMMON.VAR'
1983 c$$$      include 'COMMON.CHAIN'
1984 c$$$      include 'COMMON.INTERACT'
1985 c$$$      include 'COMMON.SBRIDGE'
1986 c$$$
1987 c$$$      double precision orig_ss_dc,orig_ss_var,orig_ss_dist
1988 c$$$      common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
1989 c$$$     +     orig_ss_dist(maxres2,maxres2)
1990 c$$$
1991 c$$$      integer n
1992 c$$$      double precision x(maxres6)
1993 c$$$      integer nf
1994 c$$$      double precision f
1995 c$$$      integer uiparm(2)
1996 c$$$      real*8 urparm(1)
1997 c$$$      external ufparm
1998 c$$$      double precision ufparm
1999 c$$$
2000 c$$$      external dist
2001 c$$$      double precision dist
2002 c$$$
2003 c$$$      integer i,j,k,ss_i,ss_j
2004 c$$$      double precision tempf,var(maxvar)
2005 c$$$
2006 c$$$
2007 c$$$      ss_i=uiparm(1)
2008 c$$$      ss_j=uiparm(2)
2009 c$$$      f=0.0D0
2010 c$$$
2011 c$$$      k=0
2012 c$$$      do i=1,nres-1
2013 c$$$        do j=1,3
2014 c$$$          k=k+1
2015 c$$$          dc(j,i)=x(k)
2016 c$$$        enddo
2017 c$$$      enddo
2018 c$$$      do i=2,nres-1
2019 c$$$        if (ialph(i,1).gt.0) then
2020 c$$$        do j=1,3
2021 c$$$          k=k+1
2022 c$$$          dc(j,i+nres)=x(k)
2023 c$$$        enddo
2024 c$$$        endif
2025 c$$$      enddo
2026 c$$$      call chainbuild_cart
2027 c$$$
2028 c$$$      call geom_to_var(nvar,var)
2029 c$$$
2030 c$$$c     Constraints on all angles
2031 c$$$      do i=1,nvar
2032 c$$$        tempf=var(i)-orig_ss_var(i)
2033 c$$$        f=f+tempf*tempf
2034 c$$$      enddo
2035 c$$$
2036 c$$$c     Constraints on all distances
2037 c$$$      do i=1,nres-1
2038 c$$$        if (i.gt.1) then
2039 c$$$          tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
2040 c$$$          f=f+tempf*tempf
2041 c$$$        endif
2042 c$$$        do j=i+1,nres
2043 c$$$          tempf=dist(j,i)-orig_ss_dist(j,i)
2044 c$$$          if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
2045 c$$$          tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
2046 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
2047 c$$$          tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
2048 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
2049 c$$$          tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
2050 c$$$          if (tempf.lt.0.0D0) f=f+tempf*tempf
2051 c$$$        enddo
2052 c$$$      enddo
2053 c$$$
2054 c$$$c     Constraints for the relevant CYS-CYS
2055 c$$$      tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
2056 c$$$      f=f+tempf*tempf
2057 c$$$CCCCCCCCCCCCCCCCC      ADD SOME ANGULAR STUFF
2058 c$$$
2059 c$$$c$$$      if (nf.ne.nfl) then
2060 c$$$c$$$        write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
2061 c$$$c$$$     +       f,dist(5+nres,14+nres)
2062 c$$$c$$$      endif
2063 c$$$
2064 c$$$      nfl=nf
2065 c$$$
2066 c$$$      return                                                            
2067 c$$$      end                                                               
2068 c$$$
2069 c$$$C-----------------------------------------------------------------------------