4c4 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 5a6 > #ifndef ISNAN 6a8 > #endif 83a86,89 > C > C 21/5/07 Calculate local sicdechain correlation energy > C > call eback_sc_corr(esccor) 99,125d104 < C call multibody(ecorr) < C < C Sum the energies < C < C scale large componenets < c#ifdef SCALE < c ecorr5_scal=1000.0 < c eel_loc_scal=100.0 < c eello_turn3_scal=100.0 < c eello_turn4_scal=100.0 < c eturn6_scal=1000.0 < c ecorr6_scal=1000.0 < c#else < c ecorr5_scal=1.0 < c eel_loc_scal=1.0 < c eello_turn3_scal=1.0 < c eello_turn4_scal=1.0 < c eturn6_scal=1.0 < c ecorr6_scal=1.0 < c#endif < c < c ecorr5=ecorr5/ecorr5_scal < c eel_loc=eel_loc/eel_loc_scal < c eello_turn3=eello_turn3/eello_turn3_scal < c eello_turn4=eello_turn4/eello_turn4_scal < c eturn6=eturn6/eturn6_scal < c ecorr6=ecorr6/ecorr6_scal 133c112 < & +wbond*estr --- > & +wbond*estr+wsccor*fact(1)*esccor 141c120 < & +wbond*estr --- > & +wbond*estr+wsccor*fact(1)*esccor 172c151,152 < energia(19)=edihcnstr --- > energia(19)=esccor > energia(20)=edihcnstr 173a154,160 > #ifdef ISNAN > #ifdef AIX > if (isnan(etot).ne.0) energia(0)=1.0d+99 > #else > if (isnan(etot)) energia(0)=1.0d+99 > #endif > #else 180a168 > #endif 201c189,190 < & wturn6*fact(5)*gcorr6_turn(j,i) --- > & wturn6*fact(5)*gcorr6_turn(j,i)+ > & wsccor*fact(2)*gsccorc(j,i) 204c193,194 < & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i) --- > & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ > & wsccor*fact(2)*gsccorx(j,i) 218c208,209 < & wturn6*fact(5)*gcorr6_turn(j,i) --- > & wturn6*fact(5)*gcorr6_turn(j,i)+ > & wsccor*fact(2)*gsccorc(j,i) 221c212,213 < & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i) --- > & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ > & wsccor*fact(1)*gsccorx(j,i) 224,225d215 < cd print '(i3,9(1pe12.4))',i,(gvdwc(k,i),k=1,3),(gelc(k,i),k=1,3), < cd & (gradc(k,i),k=1,3) 230d219 < cd write (iout,*) i,g_corr5_loc(i) 237a227 > & +wsccor*fact(1)*gsccor_loc(i) 240,244d229 < cd print*,evdw,wsc,evdw2,wscp,ees+evdw1,welec,ebe,wang, < cd & escloc,wscloc,etors,wtor,ehpb,wstrain,nss,ebr,etot < cd call enerprint(energia(0),fact) < cd call intout < cd stop 251c236 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 279c264,265 < edihcnstr=energia(19) --- > esccor=energia(19) > edihcnstr=energia(20) 289c275 < & edihcnstr,ebr*nss,etot --- > & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot 308a295 > & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ 318c305,306 < & eello_turn6,wturn6*fact(5),edihcnstr,ebr*nss,etot --- > & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor, > & edihcnstr,ebr*nss,etot 336a325 > & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ 351c340 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 360a350 > include 'COMMON.ENEPS' 368a359,363 > do i=1,210 > do j=1,2 > eneps_temp(j,i)=0.0d0 > enddo > enddo 398a394,395 > eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij) > eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij 512c509 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 519a517 > include 'COMMON.ENEPS' 526a525,529 > do i=1,210 > do j=1,2 > eneps_temp(j,i)=0.0d0 > enddo > enddo 553a557,559 > eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm) > & /dabs(eps(itypi,itypj)) > eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps(itypi,itypj) 601c607 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 609a616 > include 'COMMON.ENEPS' 616a624,628 > do i=1,210 > do j=1,2 > eneps_temp(j,i)=0.0d0 > enddo > enddo 688a701,703 > eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux > & /dabs(eps(itypi,itypj)) > eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj) 728c743 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 736a752 > include 'COMMON.ENEPS' 742a759,763 > do i=1,210 > do j=1,2 > eneps_temp(j,i)=0.0d0 > enddo > enddo 819a841,843 > eneps_temp(1,ij)=eneps_temp(1,ij)+aux*e1 > & /dabs(eps(itypi,itypj)) > eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj) 859c883 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 867a892 > include 'COMMON.ENEPS' 873a899,903 > do i=1,210 > do j=1,2 > eneps_temp(j,i)=0.0d0 > enddo > enddo 952a983,987 > eneps_temp(1,ij)=eneps_temp(1,ij)+aux*(e1+e_augm) > & /dabs(eps(itypi,itypj)) > eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj) > c eneps_temp(ij)=eneps_temp(ij) > c & +(evdwij+e_augm)/eps(itypi,itypj) 1035c1070 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 1073c1108 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 1232c1267 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 1415c1450 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 1500c1535 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 1683c1718 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 2432c2467 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 2699c2734 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 2810c2845 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 2887c2922 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 2968c3003 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 2978a3014 > double precision u(3),ud(3) 2988,2989c3024,3027 < estr=AKP*estr < c write (iout,*) "estr",estr --- > estr=0.5d0*AKP*estr > c > c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included > c 2993,2999c3031,3070 < diff=vbld(i+nres)-vbldsc0(iti) < c write (iout,*) i,iti,vbld(i+nres),vbldsc0(iti),diff, < c & AKSC(iti)*diff*diff < estr=estr+AKSC(iti)*diff*diff < do j=1,3 < gradbx(j,i)=AKSC(iti)*diff*dc(j,i+nres)/vbld(i+nres) < enddo --- > nbi=nbondterm(iti) > if (nbi.eq.1) then > diff=vbld(i+nres)-vbldsc0(1,iti) > c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, > c & AKSC(1,iti),AKSC(1,iti)*diff*diff > estr=estr+0.5d0*AKSC(1,iti)*diff*diff > do j=1,3 > gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) > enddo > else > do j=1,nbi > diff=vbld(i+nres)-vbldsc0(j,iti) > ud(j)=aksc(j,iti)*diff > u(j)=abond0(j,iti)+0.5d0*ud(j)*diff > enddo > uprod=u(1) > do j=2,nbi > uprod=uprod*u(j) > enddo > usum=0.0d0 > usumsqder=0.0d0 > do j=1,nbi > uprod1=1.0d0 > uprod2=1.0d0 > do k=1,nbi > if (k.ne.j) then > uprod1=uprod1*u(k) > uprod2=uprod2*u(k)*u(k) > endif > enddo > usum=usum+uprod1 > usumsqder=usumsqder+ud(j)*uprod2 > enddo > c write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti), > c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) > estr=estr+uprod/usum > do j=1,3 > gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) > enddo > endif 3002,3003d3072 < c write (iout,*) "estr",estr < estr=0.5d0*estr 3005a3075 > #ifdef CRYST_THETA 3014c3084 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 3244a3315,3509 > #else > C-------------------------------------------------------------------------- > subroutine ebend(etheta) > C > C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral > C angles gamma and its derivatives in consecutive thetas and gammas. > C ab initio-derived potentials from > c Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 > C > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > include 'DIMENSIONS.ZSCOPT' > include 'COMMON.LOCAL' > include 'COMMON.GEO' > include 'COMMON.INTERACT' > include 'COMMON.DERIV' > include 'COMMON.VAR' > include 'COMMON.CHAIN' > include 'COMMON.IOUNITS' > include 'COMMON.NAMES' > include 'COMMON.FFIELD' > include 'COMMON.CONTROL' > double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), > & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), > & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), > & sinph1ph2(maxdouble,maxdouble) > logical lprn /.false./, lprn1 /.false./ > etheta=0.0D0 > c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) > do i=ithet_start,ithet_end > dethetai=0.0d0 > dephii=0.0d0 > dephii1=0.0d0 > theti2=0.5d0*theta(i) > ityp2=ithetyp(itype(i-1)) > do k=1,nntheterm > coskt(k)=dcos(k*theti2) > sinkt(k)=dsin(k*theti2) > enddo > if (i.gt.3) then > #ifdef OSF > phii=phi(i) > if (phii.ne.phii) phii=150.0 > #else > phii=phi(i) > #endif > ityp1=ithetyp(itype(i-2)) > do k=1,nsingle > cosph1(k)=dcos(k*phii) > sinph1(k)=dsin(k*phii) > enddo > else > phii=0.0d0 > ityp1=nthetyp+1 > do k=1,nsingle > cosph1(k)=0.0d0 > sinph1(k)=0.0d0 > enddo > endif > if (i.lt.nres) then > #ifdef OSF > phii1=phi(i+1) > if (phii1.ne.phii1) phii1=150.0 > phii1=pinorm(phii1) > #else > phii1=phi(i+1) > #endif > ityp3=ithetyp(itype(i)) > do k=1,nsingle > cosph2(k)=dcos(k*phii1) > sinph2(k)=dsin(k*phii1) > enddo > else > phii1=0.0d0 > ityp3=nthetyp+1 > do k=1,nsingle > cosph2(k)=0.0d0 > sinph2(k)=0.0d0 > enddo > endif > c write (iout,*) "i",i," ityp1",itype(i-2),ityp1, > c & " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3 > c call flush(iout) > ethetai=aa0thet(ityp1,ityp2,ityp3) > do k=1,ndouble > do l=1,k-1 > ccl=cosph1(l)*cosph2(k-l) > ssl=sinph1(l)*sinph2(k-l) > scl=sinph1(l)*cosph2(k-l) > csl=cosph1(l)*sinph2(k-l) > cosph1ph2(l,k)=ccl-ssl > cosph1ph2(k,l)=ccl+ssl > sinph1ph2(l,k)=scl+csl > sinph1ph2(k,l)=scl-csl > enddo > enddo > if (lprn) then > write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2, > & " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1 > write (iout,*) "coskt and sinkt" > do k=1,nntheterm > write (iout,*) k,coskt(k),sinkt(k) > enddo > endif > do k=1,ntheterm > ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k) > dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3) > & *coskt(k) > if (lprn) > & write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3), > & " ethetai",ethetai > enddo > if (lprn) then > write (iout,*) "cosph and sinph" > do k=1,nsingle > write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k) > enddo > write (iout,*) "cosph1ph2 and sinph2ph2" > do k=2,ndouble > do l=1,k-1 > write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l), > & sinph1ph2(l,k),sinph1ph2(k,l) > enddo > enddo > write(iout,*) "ethetai",ethetai > endif > do m=1,ntheterm2 > do k=1,nsingle > aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k) > & +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k) > & +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k) > & +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k) > ethetai=ethetai+sinkt(m)*aux > dethetai=dethetai+0.5d0*m*aux*coskt(m) > dephii=dephii+k*sinkt(m)*( > & ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)- > & bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)) > dephii1=dephii1+k*sinkt(m)*( > & eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)- > & ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k)) > if (lprn) > & write (iout,*) "m",m," k",k," bbthet", > & bbthet(k,m,ityp1,ityp2,ityp3)," ccthet", > & ccthet(k,m,ityp1,ityp2,ityp3)," ddthet", > & ddthet(k,m,ityp1,ityp2,ityp3)," eethet", > & eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai > enddo > enddo > if (lprn) > & write(iout,*) "ethetai",ethetai > do m=1,ntheterm3 > do k=2,ndouble > do l=1,k-1 > aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+ > & ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+ > & ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+ > & ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l) > ethetai=ethetai+sinkt(m)*aux > dethetai=dethetai+0.5d0*m*coskt(m)*aux > dephii=dephii+l*sinkt(m)*( > & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)- > & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+ > & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+ > & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) > dephii1=dephii1+(k-l)*sinkt(m)*( > & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+ > & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+ > & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)- > & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) > if (lprn) then > write (iout,*) "m",m," k",k," l",l," ffthet", > & ffthet(l,k,m,ityp1,ityp2,ityp3), > & ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet", > & ggthet(l,k,m,ityp1,ityp2,ityp3), > & ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai > write (iout,*) cosph1ph2(l,k)*sinkt(m), > & cosph1ph2(k,l)*sinkt(m), > & sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m) > endif > enddo > enddo > enddo > 10 continue > if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') > & i,theta(i)*rad2deg,phii*rad2deg, > & phii1*rad2deg,ethetai > etheta=etheta+ethetai > if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii > if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 > gloc(nphi+i-2,icg)=wang*dethetai > enddo > return > end > #endif > #ifdef CRYST_SC 3252c3517 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 3525a3791,4113 > #else > c---------------------------------------------------------------------------------- > subroutine esc(escloc) > C Calculate the local energy of a side chain and its derivatives in the > C corresponding virtual-bond valence angles THETA and the spherical angles > C ALPHA and OMEGA derived from AM1 all-atom calculations. > C added by Urszula Kozlowska. 07/11/2007 > C > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > include 'DIMENSIONS.ZSCOPT' > include 'COMMON.GEO' > include 'COMMON.LOCAL' > include 'COMMON.VAR' > include 'COMMON.SCROT' > include 'COMMON.INTERACT' > include 'COMMON.DERIV' > include 'COMMON.CHAIN' > include 'COMMON.IOUNITS' > include 'COMMON.NAMES' > include 'COMMON.FFIELD' > include 'COMMON.CONTROL' > include 'COMMON.VECTORS' > double precision x_prime(3),y_prime(3),z_prime(3) > & , sumene,dsc_i,dp2_i,x(65), > & xx,yy,zz,sumene1,sumene2,sumene3,sumene4,s1,s1_6,s2,s2_6, > & de_dxx,de_dyy,de_dzz,de_dt > double precision s1_t,s1_6_t,s2_t,s2_6_t > double precision > & dXX_Ci1(3),dYY_Ci1(3),dZZ_Ci1(3),dXX_Ci(3), > & dYY_Ci(3),dZZ_Ci(3),dXX_XYZ(3),dYY_XYZ(3),dZZ_XYZ(3), > & dt_dCi(3),dt_dCi1(3) > common /sccalc/ time11,time12,time112,theti,it,nlobit > delta=0.02d0*pi > escloc=0.0D0 > do i=loc_start,loc_end > costtab(i+1) =dcos(theta(i+1)) > sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1)) > cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1))) > sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1))) > cosfac2=0.5d0/(1.0d0+costtab(i+1)) > cosfac=dsqrt(cosfac2) > sinfac2=0.5d0/(1.0d0-costtab(i+1)) > sinfac=dsqrt(sinfac2) > it=itype(i) > if (it.eq.10) goto 1 > c > C Compute the axes of tghe local cartesian coordinates system; store in > c x_prime, y_prime and z_prime > c > do j=1,3 > x_prime(j) = 0.00 > y_prime(j) = 0.00 > z_prime(j) = 0.00 > enddo > C write(2,*) "dc_norm", dc_norm(1,i+nres),dc_norm(2,i+nres), > C & dc_norm(3,i+nres) > do j = 1,3 > x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac > y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac > enddo > do j = 1,3 > z_prime(j) = -uz(j,i-1) > enddo > c write (2,*) "i",i > c write (2,*) "x_prime",(x_prime(j),j=1,3) > c write (2,*) "y_prime",(y_prime(j),j=1,3) > c write (2,*) "z_prime",(z_prime(j),j=1,3) > c write (2,*) "xx",scalar(x_prime(1),x_prime(1)), > c & " xy",scalar(x_prime(1),y_prime(1)), > c & " xz",scalar(x_prime(1),z_prime(1)), > c & " yy",scalar(y_prime(1),y_prime(1)), > c & " yz",scalar(y_prime(1),z_prime(1)), > c & " zz",scalar(z_prime(1),z_prime(1)) > c > C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i), > C to local coordinate system. Store in xx, yy, zz. > c > xx=0.0d0 > yy=0.0d0 > zz=0.0d0 > do j = 1,3 > xx = xx + x_prime(j)*dc_norm(j,i+nres) > yy = yy + y_prime(j)*dc_norm(j,i+nres) > zz = zz + z_prime(j)*dc_norm(j,i+nres) > enddo > > xxtab(i)=xx > yytab(i)=yy > zztab(i)=zz > C > C Compute the energy of the ith side cbain > C > c write (2,*) "xx",xx," yy",yy," zz",zz > it=itype(i) > do j = 1,65 > x(j) = sc_parmin(j,it) > enddo > #ifdef CHECK_COORD > Cc diagnostics - remove later > xx1 = dcos(alph(2)) > yy1 = dsin(alph(2))*dcos(omeg(2)) > zz1 = -dsin(alph(2))*dsin(omeg(2)) > write(2,'(3f8.1,3f9.3,1x,3f9.3)') > & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz, > & xx1,yy1,zz1 > C," --- ", xx_w,yy_w,zz_w > c end diagnostics > #endif > sumene1= x(1)+ x(2)*xx+ x(3)*yy+ x(4)*zz+ x(5)*xx**2 > & + x(6)*yy**2+ x(7)*zz**2+ x(8)*xx*zz+ x(9)*xx*yy > & + x(10)*yy*zz > sumene2= x(11) + x(12)*xx + x(13)*yy + x(14)*zz + x(15)*xx**2 > & + x(16)*yy**2 + x(17)*zz**2 + x(18)*xx*zz + x(19)*xx*yy > & + x(20)*yy*zz > sumene3= x(21) +x(22)*xx +x(23)*yy +x(24)*zz +x(25)*xx**2 > & +x(26)*yy**2 +x(27)*zz**2 +x(28)*xx*zz +x(29)*xx*yy > & +x(30)*yy*zz +x(31)*xx**3 +x(32)*yy**3 +x(33)*zz**3 > & +x(34)*(xx**2)*yy +x(35)*(xx**2)*zz +x(36)*(yy**2)*xx > & +x(37)*(yy**2)*zz +x(38)*(zz**2)*xx +x(39)*(zz**2)*yy > & +x(40)*xx*yy*zz > sumene4= x(41) +x(42)*xx +x(43)*yy +x(44)*zz +x(45)*xx**2 > & +x(46)*yy**2 +x(47)*zz**2 +x(48)*xx*zz +x(49)*xx*yy > & +x(50)*yy*zz +x(51)*xx**3 +x(52)*yy**3 +x(53)*zz**3 > & +x(54)*(xx**2)*yy +x(55)*(xx**2)*zz +x(56)*(yy**2)*xx > & +x(57)*(yy**2)*zz +x(58)*(zz**2)*xx +x(59)*(zz**2)*yy > & +x(60)*xx*yy*zz > dsc_i = 0.743d0+x(61) > dp2_i = 1.9d0+x(62) > dscp1=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i > & *(xx*cost2tab(i+1)+yy*sint2tab(i+1))) > dscp2=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i > & *(xx*cost2tab(i+1)-yy*sint2tab(i+1))) > s1=(1+x(63))/(0.1d0 + dscp1) > s1_6=(1+x(64))/(0.1d0 + dscp1**6) > s2=(1+x(65))/(0.1d0 + dscp2) > s2_6=(1+x(65))/(0.1d0 + dscp2**6) > sumene = ( sumene3*sint2tab(i+1) + sumene1)*(s1+s1_6) > & + (sumene4*cost2tab(i+1) +sumene2)*(s2+s2_6) > c write(2,'(i2," sumene",7f9.3)') i,sumene1,sumene2,sumene3, > c & sumene4, > c & dscp1,dscp2,sumene > c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) > escloc = escloc + sumene > c write (2,*) "escloc",escloc > if (.not. calc_grad) goto 1 > #ifdef DEBUG > C > C This section to check the numerical derivatives of the energy of ith side > C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert > C #define DEBUG in the code to turn it on. > C > write (2,*) "sumene =",sumene > aincr=1.0d-7 > xxsave=xx > xx=xx+aincr > write (2,*) xx,yy,zz > sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) > de_dxx_num=(sumenep-sumene)/aincr > xx=xxsave > write (2,*) "xx+ sumene from enesc=",sumenep > yysave=yy > yy=yy+aincr > write (2,*) xx,yy,zz > sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) > de_dyy_num=(sumenep-sumene)/aincr > yy=yysave > write (2,*) "yy+ sumene from enesc=",sumenep > zzsave=zz > zz=zz+aincr > write (2,*) xx,yy,zz > sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) > de_dzz_num=(sumenep-sumene)/aincr > zz=zzsave > write (2,*) "zz+ sumene from enesc=",sumenep > costsave=cost2tab(i+1) > sintsave=sint2tab(i+1) > cost2tab(i+1)=dcos(0.5d0*(theta(i+1)+aincr)) > sint2tab(i+1)=dsin(0.5d0*(theta(i+1)+aincr)) > sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) > de_dt_num=(sumenep-sumene)/aincr > write (2,*) " t+ sumene from enesc=",sumenep > cost2tab(i+1)=costsave > sint2tab(i+1)=sintsave > C End of diagnostics section. > #endif > C > C Compute the gradient of esc > C > pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2 > pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2 > pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2 > pom_s26=6*(1.0d0+x(65))/(0.1d0 + dscp2**6)**2 > pom_dx=dsc_i*dp2_i*cost2tab(i+1) > pom_dy=dsc_i*dp2_i*sint2tab(i+1) > pom_dt1=-0.5d0*dsc_i*dp2_i*(xx*sint2tab(i+1)-yy*cost2tab(i+1)) > pom_dt2=-0.5d0*dsc_i*dp2_i*(xx*sint2tab(i+1)+yy*cost2tab(i+1)) > pom1=(sumene3*sint2tab(i+1)+sumene1) > & *(pom_s1/dscp1+pom_s16*dscp1**4) > pom2=(sumene4*cost2tab(i+1)+sumene2) > & *(pom_s2/dscp2+pom_s26*dscp2**4) > sumene1x=x(2)+2*x(5)*xx+x(8)*zz+ x(9)*yy > sumene3x=x(22)+2*x(25)*xx+x(28)*zz+x(29)*yy+3*x(31)*xx**2 > & +2*x(34)*xx*yy +2*x(35)*xx*zz +x(36)*(yy**2) +x(38)*(zz**2) > & +x(40)*yy*zz > sumene2x=x(12)+2*x(15)*xx+x(18)*zz+ x(19)*yy > sumene4x=x(42)+2*x(45)*xx +x(48)*zz +x(49)*yy +3*x(51)*xx**2 > & +2*x(54)*xx*yy+2*x(55)*xx*zz+x(56)*(yy**2)+x(58)*(zz**2) > & +x(60)*yy*zz > de_dxx =(sumene1x+sumene3x*sint2tab(i+1))*(s1+s1_6) > & +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6) > & +(pom1+pom2)*pom_dx > #ifdef DEBUG > write(2,*), "de_dxx = ", de_dxx,de_dxx_num > #endif > C > sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz > sumene3y=x(23) +2*x(26)*yy +x(29)*xx +x(30)*zz +3*x(32)*yy**2 > & +x(34)*(xx**2) +2*x(36)*yy*xx +2*x(37)*yy*zz +x(39)*(zz**2) > & +x(40)*xx*zz > sumene2y=x(13) + 2*x(16)*yy + x(19)*xx + x(20)*zz > sumene4y=x(43)+2*x(46)*yy+x(49)*xx +x(50)*zz > & +3*x(52)*yy**2+x(54)*xx**2+2*x(56)*yy*xx +2*x(57)*yy*zz > & +x(59)*zz**2 +x(60)*xx*zz > de_dyy =(sumene1y+sumene3y*sint2tab(i+1))*(s1+s1_6) > & +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6) > & +(pom1-pom2)*pom_dy > #ifdef DEBUG > write(2,*), "de_dyy = ", de_dyy,de_dyy_num > #endif > C > de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy > & +3*x(33)*zz**2 +x(35)*xx**2 +x(37)*yy**2 +2*x(38)*zz*xx > & +2*x(39)*zz*yy +x(40)*xx*yy)*sint2tab(i+1)*(s1+s1_6) > & +(x(4) + 2*x(7)*zz+ x(8)*xx + x(10)*yy)*(s1+s1_6) > & +(x(44)+2*x(47)*zz +x(48)*xx +x(50)*yy +3*x(53)*zz**2 > & +x(55)*xx**2 +x(57)*(yy**2)+2*x(58)*zz*xx +2*x(59)*zz*yy > & +x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6) > & + ( x(14) + 2*x(17)*zz+ x(18)*xx + x(20)*yy)*(s2+s2_6) > #ifdef DEBUG > write(2,*), "de_dzz = ", de_dzz,de_dzz_num > #endif > C > de_dt = 0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6) > & -0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6) > & +pom1*pom_dt1+pom2*pom_dt2 > #ifdef DEBUG > write(2,*), "de_dt = ", de_dt,de_dt_num > #endif > c > C > cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres)) > cossc1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres)) > cosfac2xx=cosfac2*xx > sinfac2yy=sinfac2*yy > do k = 1,3 > dt_dCi(k) = -(dc_norm(k,i-1)+costtab(i+1)*dc_norm(k,i))* > & vbld_inv(i+1) > dt_dCi1(k)= -(dc_norm(k,i)+costtab(i+1)*dc_norm(k,i-1))* > & vbld_inv(i) > pom=(dC_norm(k,i+nres)-cossc*dC_norm(k,i))*vbld_inv(i+1) > pom1=(dC_norm(k,i+nres)-cossc1*dC_norm(k,i-1))*vbld_inv(i) > c write (iout,*) "i",i," k",k," pom",pom," pom1",pom1, > c & " dt_dCi",dt_dCi(k)," dt_dCi1",dt_dCi1(k) > c write (iout,*) "dC_norm",(dC_norm(j,i),j=1,3), > c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i) > dXX_Ci(k)=pom*cosfac-dt_dCi(k)*cosfac2xx > dXX_Ci1(k)=-pom1*cosfac-dt_dCi1(k)*cosfac2xx > dYY_Ci(k)=pom*sinfac+dt_dCi(k)*sinfac2yy > dYY_Ci1(k)=pom1*sinfac+dt_dCi1(k)*sinfac2yy > dZZ_Ci1(k)=0.0d0 > dZZ_Ci(k)=0.0d0 > do j=1,3 > dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres) > dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres) > enddo > > dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres)) > dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres)) > dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres)) > c > dt_dCi(k) = -dt_dCi(k)/sinttab(i+1) > dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1) > enddo > > do k=1,3 > dXX_Ctab(k,i)=dXX_Ci(k) > dXX_C1tab(k,i)=dXX_Ci1(k) > dYY_Ctab(k,i)=dYY_Ci(k) > dYY_C1tab(k,i)=dYY_Ci1(k) > dZZ_Ctab(k,i)=dZZ_Ci(k) > dZZ_C1tab(k,i)=dZZ_Ci1(k) > dXX_XYZtab(k,i)=dXX_XYZ(k) > dYY_XYZtab(k,i)=dYY_XYZ(k) > dZZ_XYZtab(k,i)=dZZ_XYZ(k) > enddo > > do k = 1,3 > c write (iout,*) "k",k," dxx_ci1",dxx_ci1(k)," dyy_ci1", > c & dyy_ci1(k)," dzz_ci1",dzz_ci1(k) > c write (iout,*) "k",k," dxx_ci",dxx_ci(k)," dyy_ci", > c & dyy_ci(k)," dzz_ci",dzz_ci(k) > c write (iout,*) "k",k," dt_dci",dt_dci(k)," dt_dci", > c & dt_dci(k) > c write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ", > c & dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k) > gscloc(k,i-1)=gscloc(k,i-1)+de_dxx*dxx_ci1(k) > & +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k) > gscloc(k,i)=gscloc(k,i)+de_dxx*dxx_Ci(k) > & +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k) > gsclocx(k,i)= de_dxx*dxx_XYZ(k) > & +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k) > enddo > c write(iout,*) "ENERGY GRAD = ", (gscloc(k,i-1),k=1,3), > c & (gscloc(k,i),k=1,3),(gsclocx(k,i),k=1,3) > > C to check gradient call subroutine check_grad > > 1 continue > enddo > return > end > #endif 3563c4151 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 3611c4199 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 3694c4282 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 3780c4368 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 3847a4436,4486 > subroutine eback_sc_corr(esccor) > c 7/21/2007 Correlations between the backbone-local and side-chain-local > c conformational states; temporarily implemented as differences > c between UNRES torsional potentials (dependent on three types of > c residues) and the torsional potentials dependent on all 20 types > c of residues computed from AM1 energy surfaces of terminally-blocked > c amino-acid residues. > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > include 'DIMENSIONS.ZSCOPT' > include 'COMMON.VAR' > include 'COMMON.GEO' > include 'COMMON.LOCAL' > include 'COMMON.TORSION' > include 'COMMON.SCCOR' > include 'COMMON.INTERACT' > include 'COMMON.DERIV' > include 'COMMON.CHAIN' > include 'COMMON.NAMES' > include 'COMMON.IOUNITS' > include 'COMMON.FFIELD' > include 'COMMON.CONTROL' > logical lprn > C Set lprn=.true. for debugging > lprn=.false. > c lprn=.true. > c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor > esccor=0.0D0 > do i=iphi_start,iphi_end > esccor_ii=0.0D0 > itori=itype(i-2) > itori1=itype(i-1) > phii=phi(i) > gloci=0.0D0 > do j=1,nterm_sccor > v1ij=v1sccor(j,itori,itori1) > v2ij=v2sccor(j,itori,itori1) > cosphi=dcos(j*phii) > sinphi=dsin(j*phii) > esccor=esccor+v1ij*cosphi+v2ij*sinphi > gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) > enddo > if (lprn) > & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') > & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, > & (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6) > gsccor_loc(i-3)=gloci > enddo > return > end > c------------------------------------------------------------------------------ 4003c4642 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 4189c4828 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 4498c5137 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 4565c5204 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 4942c5581 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 5059c5698 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 5460c6099 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 5597c6236 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 5703c6342 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 5890c6529 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 6006c6645 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT' 6252c6891 < include 'sizesclu.dat' --- > include 'DIMENSIONS.ZSCOPT'