2 < include 'sizesclu.dat'
4 > include 'DIMENSIONS.ZSCOPT'
11 > C 21/5/07 Calculate local sicdechain correlation energy
13 > call eback_sc_corr(esccor)
15 < C call multibody(ecorr)
19 < C scale large componenets
21 < c ecorr5_scal=1000.0
22 < c eel_loc_scal=100.0
23 < c eello_turn3_scal=100.0
24 < c eello_turn4_scal=100.0
25 < c eturn6_scal=1000.0
26 < c ecorr6_scal=1000.0
30 < c eello_turn3_scal=1.0
31 < c eello_turn4_scal=1.0
36 < c ecorr5=ecorr5/ecorr5_scal
37 < c eel_loc=eel_loc/eel_loc_scal
38 < c eello_turn3=eello_turn3/eello_turn3_scal
39 < c eello_turn4=eello_turn4/eello_turn4_scal
40 < c eturn6=eturn6/eturn6_scal
41 < c ecorr6=ecorr6/ecorr6_scal
45 > & +wbond*estr+wsccor*fact(1)*esccor
49 > & +wbond*estr+wsccor*fact(1)*esccor
51 < energia(19)=edihcnstr
54 > energia(20)=edihcnstr
58 > if (isnan(etot).ne.0) energia(0)=1.0d+99
60 > if (isnan(etot)) energia(0)=1.0d+99
66 < & wturn6*fact(5)*gcorr6_turn(j,i)
68 > & wturn6*fact(5)*gcorr6_turn(j,i)+
69 > & wsccor*fact(2)*gsccorc(j,i)
71 < & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)
73 > & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
74 > & wsccor*fact(2)*gsccorx(j,i)
76 < & wturn6*fact(5)*gcorr6_turn(j,i)
78 > & wturn6*fact(5)*gcorr6_turn(j,i)+
79 > & wsccor*fact(2)*gsccorc(j,i)
81 < & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)
83 > & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
84 > & wsccor*fact(1)*gsccorx(j,i)
86 < cd print '(i3,9(1pe12.4))',i,(gvdwc(k,i),k=1,3),(gelc(k,i),k=1,3),
87 < cd & (gradc(k,i),k=1,3)
89 < cd write (iout,*) i,g_corr5_loc(i)
91 > & +wsccor*fact(1)*gsccor_loc(i)
93 < cd print*,evdw,wsc,evdw2,wscp,ees+evdw1,welec,ebe,wang,
94 < cd & escloc,wscloc,etors,wtor,ehpb,wstrain,nss,ebr,etot
95 < cd call enerprint(energia(0),fact)
99 < include 'sizesclu.dat'
101 > include 'DIMENSIONS.ZSCOPT'
103 < edihcnstr=energia(19)
106 > edihcnstr=energia(20)
108 < & edihcnstr,ebr*nss,etot
110 > & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot
112 > & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
114 < & eello_turn6,wturn6*fact(5),edihcnstr,ebr*nss,etot
116 > & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
117 > & edihcnstr,ebr*nss,etot
119 > & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
121 < include 'sizesclu.dat'
123 > include 'DIMENSIONS.ZSCOPT'
125 > include 'COMMON.ENEPS'
129 > eneps_temp(j,i)=0.0d0
133 > eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij)
134 > eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij
136 < include 'sizesclu.dat'
138 > include 'DIMENSIONS.ZSCOPT'
140 > include 'COMMON.ENEPS'
144 > eneps_temp(j,i)=0.0d0
148 > eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm)
149 > & /dabs(eps(itypi,itypj))
150 > eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps(itypi,itypj)
152 < include 'sizesclu.dat'
154 > include 'DIMENSIONS.ZSCOPT'
156 > include 'COMMON.ENEPS'
160 > eneps_temp(j,i)=0.0d0
164 > eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux
165 > & /dabs(eps(itypi,itypj))
166 > eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj)
168 < include 'sizesclu.dat'
170 > include 'DIMENSIONS.ZSCOPT'
172 > include 'COMMON.ENEPS'
176 > eneps_temp(j,i)=0.0d0
180 > eneps_temp(1,ij)=eneps_temp(1,ij)+aux*e1
181 > & /dabs(eps(itypi,itypj))
182 > eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
184 < include 'sizesclu.dat'
186 > include 'DIMENSIONS.ZSCOPT'
188 > include 'COMMON.ENEPS'
192 > eneps_temp(j,i)=0.0d0
196 > eneps_temp(1,ij)=eneps_temp(1,ij)+aux*(e1+e_augm)
197 > & /dabs(eps(itypi,itypj))
198 > eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
199 > c eneps_temp(ij)=eneps_temp(ij)
200 > c & +(evdwij+e_augm)/eps(itypi,itypj)
202 < include 'sizesclu.dat'
204 > include 'DIMENSIONS.ZSCOPT'
206 < include 'sizesclu.dat'
208 > include 'DIMENSIONS.ZSCOPT'
210 < include 'sizesclu.dat'
212 > include 'DIMENSIONS.ZSCOPT'
214 < include 'sizesclu.dat'
216 > include 'DIMENSIONS.ZSCOPT'
218 < include 'sizesclu.dat'
220 > include 'DIMENSIONS.ZSCOPT'
222 < include 'sizesclu.dat'
224 > include 'DIMENSIONS.ZSCOPT'
226 < include 'sizesclu.dat'
228 > include 'DIMENSIONS.ZSCOPT'
230 < include 'sizesclu.dat'
232 > include 'DIMENSIONS.ZSCOPT'
234 < include 'sizesclu.dat'
236 > include 'DIMENSIONS.ZSCOPT'
238 < include 'sizesclu.dat'
240 > include 'DIMENSIONS.ZSCOPT'
242 < include 'sizesclu.dat'
244 > include 'DIMENSIONS.ZSCOPT'
246 > double precision u(3),ud(3)
249 < c write (iout,*) "estr",estr
251 > estr=0.5d0*AKP*estr
253 > c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
256 < diff=vbld(i+nres)-vbldsc0(iti)
257 < c write (iout,*) i,iti,vbld(i+nres),vbldsc0(iti),diff,
258 < c & AKSC(iti)*diff*diff
259 < estr=estr+AKSC(iti)*diff*diff
261 < gradbx(j,i)=AKSC(iti)*diff*dc(j,i+nres)/vbld(i+nres)
266 > diff=vbld(i+nres)-vbldsc0(1,iti)
267 > c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
268 > c & AKSC(1,iti),AKSC(1,iti)*diff*diff
269 > estr=estr+0.5d0*AKSC(1,iti)*diff*diff
271 > gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
275 > diff=vbld(i+nres)-vbldsc0(j,iti)
276 > ud(j)=aksc(j,iti)*diff
277 > u(j)=abond0(j,iti)+0.5d0*ud(j)*diff
291 > uprod2=uprod2*u(k)*u(k)
295 > usumsqder=usumsqder+ud(j)*uprod2
297 > c write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
298 > c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
299 > estr=estr+uprod/usum
301 > gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
305 < c write (iout,*) "estr",estr
310 < include 'sizesclu.dat'
312 > include 'DIMENSIONS.ZSCOPT'
315 > C--------------------------------------------------------------------------
316 > subroutine ebend(etheta)
318 > C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
319 > C angles gamma and its derivatives in consecutive thetas and gammas.
320 > C ab initio-derived potentials from
321 > c Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
323 > implicit real*8 (a-h,o-z)
324 > include 'DIMENSIONS'
325 > include 'DIMENSIONS.ZSCOPT'
326 > include 'COMMON.LOCAL'
327 > include 'COMMON.GEO'
328 > include 'COMMON.INTERACT'
329 > include 'COMMON.DERIV'
330 > include 'COMMON.VAR'
331 > include 'COMMON.CHAIN'
332 > include 'COMMON.IOUNITS'
333 > include 'COMMON.NAMES'
334 > include 'COMMON.FFIELD'
335 > include 'COMMON.CONTROL'
336 > double precision coskt(mmaxtheterm),sinkt(mmaxtheterm),
337 > & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle),
338 > & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble),
339 > & sinph1ph2(maxdouble,maxdouble)
340 > logical lprn /.false./, lprn1 /.false./
342 > c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
343 > do i=ithet_start,ithet_end
347 > theti2=0.5d0*theta(i)
348 > ityp2=ithetyp(itype(i-1))
350 > coskt(k)=dcos(k*theti2)
351 > sinkt(k)=dsin(k*theti2)
356 > if (phii.ne.phii) phii=150.0
360 > ityp1=ithetyp(itype(i-2))
362 > cosph1(k)=dcos(k*phii)
363 > sinph1(k)=dsin(k*phii)
373 > if (i.lt.nres) then
376 > if (phii1.ne.phii1) phii1=150.0
377 > phii1=pinorm(phii1)
381 > ityp3=ithetyp(itype(i))
383 > cosph2(k)=dcos(k*phii1)
384 > sinph2(k)=dsin(k*phii1)
394 > c write (iout,*) "i",i," ityp1",itype(i-2),ityp1,
395 > c & " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3
397 > ethetai=aa0thet(ityp1,ityp2,ityp3)
400 > ccl=cosph1(l)*cosph2(k-l)
401 > ssl=sinph1(l)*sinph2(k-l)
402 > scl=sinph1(l)*cosph2(k-l)
403 > csl=cosph1(l)*sinph2(k-l)
404 > cosph1ph2(l,k)=ccl-ssl
405 > cosph1ph2(k,l)=ccl+ssl
406 > sinph1ph2(l,k)=scl+csl
407 > sinph1ph2(k,l)=scl-csl
411 > write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,
412 > & " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1
413 > write (iout,*) "coskt and sinkt"
415 > write (iout,*) k,coskt(k),sinkt(k)
419 > ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
420 > dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
423 > & write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
424 > & " ethetai",ethetai
427 > write (iout,*) "cosph and sinph"
429 > write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k)
431 > write (iout,*) "cosph1ph2 and sinph2ph2"
434 > write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l),
435 > & sinph1ph2(l,k),sinph1ph2(k,l)
438 > write(iout,*) "ethetai",ethetai
442 > aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
443 > & +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
444 > & +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
445 > & +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
446 > ethetai=ethetai+sinkt(m)*aux
447 > dethetai=dethetai+0.5d0*m*aux*coskt(m)
448 > dephii=dephii+k*sinkt(m)*(
449 > & ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
450 > & bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
451 > dephii1=dephii1+k*sinkt(m)*(
452 > & eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
453 > & ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
455 > & write (iout,*) "m",m," k",k," bbthet",
456 > & bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
457 > & ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
458 > & ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
459 > & eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
463 > & write(iout,*) "ethetai",ethetai
467 > aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
468 > & ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
469 > & ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
470 > & ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
471 > ethetai=ethetai+sinkt(m)*aux
472 > dethetai=dethetai+0.5d0*m*coskt(m)*aux
473 > dephii=dephii+l*sinkt(m)*(
474 > & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
475 > & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
476 > & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
477 > & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
478 > dephii1=dephii1+(k-l)*sinkt(m)*(
479 > & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
480 > & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
481 > & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
482 > & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
484 > write (iout,*) "m",m," k",k," l",l," ffthet",
485 > & ffthet(l,k,m,ityp1,ityp2,ityp3),
486 > & ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
487 > & ggthet(l,k,m,ityp1,ityp2,ityp3),
488 > & ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
489 > write (iout,*) cosph1ph2(l,k)*sinkt(m),
490 > & cosph1ph2(k,l)*sinkt(m),
491 > & sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
497 > if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
498 > & i,theta(i)*rad2deg,phii*rad2deg,
499 > & phii1*rad2deg,ethetai
500 > etheta=etheta+ethetai
501 > if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
502 > if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
503 > gloc(nphi+i-2,icg)=wang*dethetai
510 < include 'sizesclu.dat'
512 > include 'DIMENSIONS.ZSCOPT'
515 > c----------------------------------------------------------------------------------
516 > subroutine esc(escloc)
517 > C Calculate the local energy of a side chain and its derivatives in the
518 > C corresponding virtual-bond valence angles THETA and the spherical angles
519 > C ALPHA and OMEGA derived from AM1 all-atom calculations.
520 > C added by Urszula Kozlowska. 07/11/2007
522 > implicit real*8 (a-h,o-z)
523 > include 'DIMENSIONS'
524 > include 'DIMENSIONS.ZSCOPT'
525 > include 'COMMON.GEO'
526 > include 'COMMON.LOCAL'
527 > include 'COMMON.VAR'
528 > include 'COMMON.SCROT'
529 > include 'COMMON.INTERACT'
530 > include 'COMMON.DERIV'
531 > include 'COMMON.CHAIN'
532 > include 'COMMON.IOUNITS'
533 > include 'COMMON.NAMES'
534 > include 'COMMON.FFIELD'
535 > include 'COMMON.CONTROL'
536 > include 'COMMON.VECTORS'
537 > double precision x_prime(3),y_prime(3),z_prime(3)
538 > & , sumene,dsc_i,dp2_i,x(65),
539 > & xx,yy,zz,sumene1,sumene2,sumene3,sumene4,s1,s1_6,s2,s2_6,
540 > & de_dxx,de_dyy,de_dzz,de_dt
541 > double precision s1_t,s1_6_t,s2_t,s2_6_t
543 > & dXX_Ci1(3),dYY_Ci1(3),dZZ_Ci1(3),dXX_Ci(3),
544 > & dYY_Ci(3),dZZ_Ci(3),dXX_XYZ(3),dYY_XYZ(3),dZZ_XYZ(3),
545 > & dt_dCi(3),dt_dCi1(3)
546 > common /sccalc/ time11,time12,time112,theti,it,nlobit
549 > do i=loc_start,loc_end
550 > costtab(i+1) =dcos(theta(i+1))
551 > sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
552 > cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
553 > sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
554 > cosfac2=0.5d0/(1.0d0+costtab(i+1))
555 > cosfac=dsqrt(cosfac2)
556 > sinfac2=0.5d0/(1.0d0-costtab(i+1))
557 > sinfac=dsqrt(sinfac2)
559 > if (it.eq.10) goto 1
561 > C Compute the axes of tghe local cartesian coordinates system; store in
562 > c x_prime, y_prime and z_prime
569 > C write(2,*) "dc_norm", dc_norm(1,i+nres),dc_norm(2,i+nres),
570 > C & dc_norm(3,i+nres)
572 > x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
573 > y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
576 > z_prime(j) = -uz(j,i-1)
578 > c write (2,*) "i",i
579 > c write (2,*) "x_prime",(x_prime(j),j=1,3)
580 > c write (2,*) "y_prime",(y_prime(j),j=1,3)
581 > c write (2,*) "z_prime",(z_prime(j),j=1,3)
582 > c write (2,*) "xx",scalar(x_prime(1),x_prime(1)),
583 > c & " xy",scalar(x_prime(1),y_prime(1)),
584 > c & " xz",scalar(x_prime(1),z_prime(1)),
585 > c & " yy",scalar(y_prime(1),y_prime(1)),
586 > c & " yz",scalar(y_prime(1),z_prime(1)),
587 > c & " zz",scalar(z_prime(1),z_prime(1))
589 > C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
590 > C to local coordinate system. Store in xx, yy, zz.
596 > xx = xx + x_prime(j)*dc_norm(j,i+nres)
597 > yy = yy + y_prime(j)*dc_norm(j,i+nres)
598 > zz = zz + z_prime(j)*dc_norm(j,i+nres)
605 > C Compute the energy of the ith side cbain
607 > c write (2,*) "xx",xx," yy",yy," zz",zz
610 > x(j) = sc_parmin(j,it)
613 > Cc diagnostics - remove later
614 > xx1 = dcos(alph(2))
615 > yy1 = dsin(alph(2))*dcos(omeg(2))
616 > zz1 = -dsin(alph(2))*dsin(omeg(2))
617 > write(2,'(3f8.1,3f9.3,1x,3f9.3)')
618 > & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
620 > C," --- ", xx_w,yy_w,zz_w
623 > sumene1= x(1)+ x(2)*xx+ x(3)*yy+ x(4)*zz+ x(5)*xx**2
624 > & + x(6)*yy**2+ x(7)*zz**2+ x(8)*xx*zz+ x(9)*xx*yy
626 > sumene2= x(11) + x(12)*xx + x(13)*yy + x(14)*zz + x(15)*xx**2
627 > & + x(16)*yy**2 + x(17)*zz**2 + x(18)*xx*zz + x(19)*xx*yy
629 > sumene3= x(21) +x(22)*xx +x(23)*yy +x(24)*zz +x(25)*xx**2
630 > & +x(26)*yy**2 +x(27)*zz**2 +x(28)*xx*zz +x(29)*xx*yy
631 > & +x(30)*yy*zz +x(31)*xx**3 +x(32)*yy**3 +x(33)*zz**3
632 > & +x(34)*(xx**2)*yy +x(35)*(xx**2)*zz +x(36)*(yy**2)*xx
633 > & +x(37)*(yy**2)*zz +x(38)*(zz**2)*xx +x(39)*(zz**2)*yy
635 > sumene4= x(41) +x(42)*xx +x(43)*yy +x(44)*zz +x(45)*xx**2
636 > & +x(46)*yy**2 +x(47)*zz**2 +x(48)*xx*zz +x(49)*xx*yy
637 > & +x(50)*yy*zz +x(51)*xx**3 +x(52)*yy**3 +x(53)*zz**3
638 > & +x(54)*(xx**2)*yy +x(55)*(xx**2)*zz +x(56)*(yy**2)*xx
639 > & +x(57)*(yy**2)*zz +x(58)*(zz**2)*xx +x(59)*(zz**2)*yy
641 > dsc_i = 0.743d0+x(61)
642 > dp2_i = 1.9d0+x(62)
643 > dscp1=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i
644 > & *(xx*cost2tab(i+1)+yy*sint2tab(i+1)))
645 > dscp2=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i
646 > & *(xx*cost2tab(i+1)-yy*sint2tab(i+1)))
647 > s1=(1+x(63))/(0.1d0 + dscp1)
648 > s1_6=(1+x(64))/(0.1d0 + dscp1**6)
649 > s2=(1+x(65))/(0.1d0 + dscp2)
650 > s2_6=(1+x(65))/(0.1d0 + dscp2**6)
651 > sumene = ( sumene3*sint2tab(i+1) + sumene1)*(s1+s1_6)
652 > & + (sumene4*cost2tab(i+1) +sumene2)*(s2+s2_6)
653 > c write(2,'(i2," sumene",7f9.3)') i,sumene1,sumene2,sumene3,
655 > c & dscp1,dscp2,sumene
656 > c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
657 > escloc = escloc + sumene
658 > c write (2,*) "escloc",escloc
659 > if (.not. calc_grad) goto 1
662 > C This section to check the numerical derivatives of the energy of ith side
663 > C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
664 > C #define DEBUG in the code to turn it on.
666 > write (2,*) "sumene =",sumene
670 > write (2,*) xx,yy,zz
671 > sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
672 > de_dxx_num=(sumenep-sumene)/aincr
674 > write (2,*) "xx+ sumene from enesc=",sumenep
677 > write (2,*) xx,yy,zz
678 > sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
679 > de_dyy_num=(sumenep-sumene)/aincr
681 > write (2,*) "yy+ sumene from enesc=",sumenep
684 > write (2,*) xx,yy,zz
685 > sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
686 > de_dzz_num=(sumenep-sumene)/aincr
688 > write (2,*) "zz+ sumene from enesc=",sumenep
689 > costsave=cost2tab(i+1)
690 > sintsave=sint2tab(i+1)
691 > cost2tab(i+1)=dcos(0.5d0*(theta(i+1)+aincr))
692 > sint2tab(i+1)=dsin(0.5d0*(theta(i+1)+aincr))
693 > sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
694 > de_dt_num=(sumenep-sumene)/aincr
695 > write (2,*) " t+ sumene from enesc=",sumenep
696 > cost2tab(i+1)=costsave
697 > sint2tab(i+1)=sintsave
698 > C End of diagnostics section.
701 > C Compute the gradient of esc
703 > pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2
704 > pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2
705 > pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2
706 > pom_s26=6*(1.0d0+x(65))/(0.1d0 + dscp2**6)**2
707 > pom_dx=dsc_i*dp2_i*cost2tab(i+1)
708 > pom_dy=dsc_i*dp2_i*sint2tab(i+1)
709 > pom_dt1=-0.5d0*dsc_i*dp2_i*(xx*sint2tab(i+1)-yy*cost2tab(i+1))
710 > pom_dt2=-0.5d0*dsc_i*dp2_i*(xx*sint2tab(i+1)+yy*cost2tab(i+1))
711 > pom1=(sumene3*sint2tab(i+1)+sumene1)
712 > & *(pom_s1/dscp1+pom_s16*dscp1**4)
713 > pom2=(sumene4*cost2tab(i+1)+sumene2)
714 > & *(pom_s2/dscp2+pom_s26*dscp2**4)
715 > sumene1x=x(2)+2*x(5)*xx+x(8)*zz+ x(9)*yy
716 > sumene3x=x(22)+2*x(25)*xx+x(28)*zz+x(29)*yy+3*x(31)*xx**2
717 > & +2*x(34)*xx*yy +2*x(35)*xx*zz +x(36)*(yy**2) +x(38)*(zz**2)
719 > sumene2x=x(12)+2*x(15)*xx+x(18)*zz+ x(19)*yy
720 > sumene4x=x(42)+2*x(45)*xx +x(48)*zz +x(49)*yy +3*x(51)*xx**2
721 > & +2*x(54)*xx*yy+2*x(55)*xx*zz+x(56)*(yy**2)+x(58)*(zz**2)
723 > de_dxx =(sumene1x+sumene3x*sint2tab(i+1))*(s1+s1_6)
724 > & +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6)
725 > & +(pom1+pom2)*pom_dx
727 > write(2,*), "de_dxx = ", de_dxx,de_dxx_num
730 > sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz
731 > sumene3y=x(23) +2*x(26)*yy +x(29)*xx +x(30)*zz +3*x(32)*yy**2
732 > & +x(34)*(xx**2) +2*x(36)*yy*xx +2*x(37)*yy*zz +x(39)*(zz**2)
734 > sumene2y=x(13) + 2*x(16)*yy + x(19)*xx + x(20)*zz
735 > sumene4y=x(43)+2*x(46)*yy+x(49)*xx +x(50)*zz
736 > & +3*x(52)*yy**2+x(54)*xx**2+2*x(56)*yy*xx +2*x(57)*yy*zz
737 > & +x(59)*zz**2 +x(60)*xx*zz
738 > de_dyy =(sumene1y+sumene3y*sint2tab(i+1))*(s1+s1_6)
739 > & +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6)
740 > & +(pom1-pom2)*pom_dy
742 > write(2,*), "de_dyy = ", de_dyy,de_dyy_num
745 > de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy
746 > & +3*x(33)*zz**2 +x(35)*xx**2 +x(37)*yy**2 +2*x(38)*zz*xx
747 > & +2*x(39)*zz*yy +x(40)*xx*yy)*sint2tab(i+1)*(s1+s1_6)
748 > & +(x(4) + 2*x(7)*zz+ x(8)*xx + x(10)*yy)*(s1+s1_6)
749 > & +(x(44)+2*x(47)*zz +x(48)*xx +x(50)*yy +3*x(53)*zz**2
750 > & +x(55)*xx**2 +x(57)*(yy**2)+2*x(58)*zz*xx +2*x(59)*zz*yy
751 > & +x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6)
752 > & + ( x(14) + 2*x(17)*zz+ x(18)*xx + x(20)*yy)*(s2+s2_6)
754 > write(2,*), "de_dzz = ", de_dzz,de_dzz_num
757 > de_dt = 0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6)
758 > & -0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6)
759 > & +pom1*pom_dt1+pom2*pom_dt2
761 > write(2,*), "de_dt = ", de_dt,de_dt_num
765 > cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
766 > cossc1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
767 > cosfac2xx=cosfac2*xx
768 > sinfac2yy=sinfac2*yy
770 > dt_dCi(k) = -(dc_norm(k,i-1)+costtab(i+1)*dc_norm(k,i))*
772 > dt_dCi1(k)= -(dc_norm(k,i)+costtab(i+1)*dc_norm(k,i-1))*
774 > pom=(dC_norm(k,i+nres)-cossc*dC_norm(k,i))*vbld_inv(i+1)
775 > pom1=(dC_norm(k,i+nres)-cossc1*dC_norm(k,i-1))*vbld_inv(i)
776 > c write (iout,*) "i",i," k",k," pom",pom," pom1",pom1,
777 > c & " dt_dCi",dt_dCi(k)," dt_dCi1",dt_dCi1(k)
778 > c write (iout,*) "dC_norm",(dC_norm(j,i),j=1,3),
779 > c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
780 > dXX_Ci(k)=pom*cosfac-dt_dCi(k)*cosfac2xx
781 > dXX_Ci1(k)=-pom1*cosfac-dt_dCi1(k)*cosfac2xx
782 > dYY_Ci(k)=pom*sinfac+dt_dCi(k)*sinfac2yy
783 > dYY_Ci1(k)=pom1*sinfac+dt_dCi1(k)*sinfac2yy
787 > dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
788 > dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
791 > dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
792 > dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres))
793 > dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres))
795 > dt_dCi(k) = -dt_dCi(k)/sinttab(i+1)
796 > dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1)
800 > dXX_Ctab(k,i)=dXX_Ci(k)
801 > dXX_C1tab(k,i)=dXX_Ci1(k)
802 > dYY_Ctab(k,i)=dYY_Ci(k)
803 > dYY_C1tab(k,i)=dYY_Ci1(k)
804 > dZZ_Ctab(k,i)=dZZ_Ci(k)
805 > dZZ_C1tab(k,i)=dZZ_Ci1(k)
806 > dXX_XYZtab(k,i)=dXX_XYZ(k)
807 > dYY_XYZtab(k,i)=dYY_XYZ(k)
808 > dZZ_XYZtab(k,i)=dZZ_XYZ(k)
812 > c write (iout,*) "k",k," dxx_ci1",dxx_ci1(k)," dyy_ci1",
813 > c & dyy_ci1(k)," dzz_ci1",dzz_ci1(k)
814 > c write (iout,*) "k",k," dxx_ci",dxx_ci(k)," dyy_ci",
815 > c & dyy_ci(k)," dzz_ci",dzz_ci(k)
816 > c write (iout,*) "k",k," dt_dci",dt_dci(k)," dt_dci",
818 > c write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ",
819 > c & dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k)
820 > gscloc(k,i-1)=gscloc(k,i-1)+de_dxx*dxx_ci1(k)
821 > & +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k)
822 > gscloc(k,i)=gscloc(k,i)+de_dxx*dxx_Ci(k)
823 > & +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k)
824 > gsclocx(k,i)= de_dxx*dxx_XYZ(k)
825 > & +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k)
827 > c write(iout,*) "ENERGY GRAD = ", (gscloc(k,i-1),k=1,3),
828 > c & (gscloc(k,i),k=1,3),(gsclocx(k,i),k=1,3)
830 > C to check gradient call subroutine check_grad
838 < include 'sizesclu.dat'
840 > include 'DIMENSIONS.ZSCOPT'
842 < include 'sizesclu.dat'
844 > include 'DIMENSIONS.ZSCOPT'
846 < include 'sizesclu.dat'
848 > include 'DIMENSIONS.ZSCOPT'
850 < include 'sizesclu.dat'
852 > include 'DIMENSIONS.ZSCOPT'
854 > subroutine eback_sc_corr(esccor)
855 > c 7/21/2007 Correlations between the backbone-local and side-chain-local
856 > c conformational states; temporarily implemented as differences
857 > c between UNRES torsional potentials (dependent on three types of
858 > c residues) and the torsional potentials dependent on all 20 types
859 > c of residues computed from AM1 energy surfaces of terminally-blocked
860 > c amino-acid residues.
861 > implicit real*8 (a-h,o-z)
862 > include 'DIMENSIONS'
863 > include 'DIMENSIONS.ZSCOPT'
864 > include 'COMMON.VAR'
865 > include 'COMMON.GEO'
866 > include 'COMMON.LOCAL'
867 > include 'COMMON.TORSION'
868 > include 'COMMON.SCCOR'
869 > include 'COMMON.INTERACT'
870 > include 'COMMON.DERIV'
871 > include 'COMMON.CHAIN'
872 > include 'COMMON.NAMES'
873 > include 'COMMON.IOUNITS'
874 > include 'COMMON.FFIELD'
875 > include 'COMMON.CONTROL'
877 > C Set lprn=.true. for debugging
880 > c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
882 > do i=iphi_start,iphi_end
889 > v1ij=v1sccor(j,itori,itori1)
890 > v2ij=v2sccor(j,itori,itori1)
891 > cosphi=dcos(j*phii)
892 > sinphi=dsin(j*phii)
893 > esccor=esccor+v1ij*cosphi+v2ij*sinphi
894 > gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
897 > & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
898 > & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
899 > & (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
900 > gsccor_loc(i-3)=gloci
904 > c------------------------------------------------------------------------------
906 < include 'sizesclu.dat'
908 > include 'DIMENSIONS.ZSCOPT'
910 < include 'sizesclu.dat'
912 > include 'DIMENSIONS.ZSCOPT'
914 < include 'sizesclu.dat'
916 > include 'DIMENSIONS.ZSCOPT'
918 < include 'sizesclu.dat'
920 > include 'DIMENSIONS.ZSCOPT'
922 < include 'sizesclu.dat'
924 > include 'DIMENSIONS.ZSCOPT'
926 < include 'sizesclu.dat'
928 > include 'DIMENSIONS.ZSCOPT'
930 < include 'sizesclu.dat'
932 > include 'DIMENSIONS.ZSCOPT'
934 < include 'sizesclu.dat'
936 > include 'DIMENSIONS.ZSCOPT'
938 < include 'sizesclu.dat'
940 > include 'DIMENSIONS.ZSCOPT'
942 < include 'sizesclu.dat'
944 > include 'DIMENSIONS.ZSCOPT'
946 < include 'sizesclu.dat'
948 > include 'DIMENSIONS.ZSCOPT'
950 < include 'sizesclu.dat'
952 > include 'DIMENSIONS.ZSCOPT'