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