C Calculate the virtual-bond-angle energy.
C
if (wang.gt.0d0) then
+ if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
call ebend(ebe,ethetacnstr)
+ endif
+C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
+C energy function
+ if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
+ call ebend_kcc(ebe,ethetacnstr)
+ endif
else
ebe=0
ethetacnstr=0
C Calculate the virtual-bond torsional energy.
C
cd print *,'nterm=',nterm
+C print *,"tor",tor_mode
if (wtor.gt.0) then
+ if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
call etor(etors,edihcnstr)
+ endif
+C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
+C energy function
+ if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
+ call etor_kcc(etors,edihcnstr)
+ endif
else
etors=0
edihcnstr=0
C
C 6/23/01 Calculate double-torsional energy
C
- if (wtor_d.gt.0) then
+ if ((wtor_d.gt.0).and.((tor_mode.eq.0).or.(tor_mode.eq.2))) then
call etor_d(etors_d)
else
etors_d=0
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
& +welec*gshieldc(j,i)
+ & +wcorr*gshieldc_ec(j,i)
+ & +wturn3*gshieldc_t3(j,i)
+ & +wturn4*gshieldc_t4(j,i)
+ & +wel_loc*gshieldc_ll(j,i)
+
enddo
enddo
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
& +welec*gshieldc(j,i)
+ & +wcorr*gshieldc_ec(j,i)
+ & +wturn4*gshieldc_t4(j,i)
+ & +wel_loc*gshieldc_ll(j,i)
+
enddo
enddo
& +gradafm(j,i)
& +welec*gshieldc(j,i)
& +welec*gshieldc_loc(j,i)
+ & +wcorr*gshieldc_ec(j,i)
+ & +wcorr*gshieldc_loc_ec(j,i)
+ & +wturn3*gshieldc_t3(j,i)
+ & +wturn3*gshieldc_loc_t3(j,i)
+ & +wturn4*gshieldc_t4(j,i)
+ & +wturn4*gshieldc_loc_t4(j,i)
+ & +wel_loc*gshieldc_ll(j,i)
+ & +wel_loc*gshieldc_loc_ll(j,i)
+
+
+
+
#else
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& 0.5d0*(wscp*gvdwc_scpp(j,i)+
- & welec*gelc_long(j,i)
+ & welec*gelc_long(j,i)+
& wel_loc*gel_loc_long(j,i)+
& wcorr*gcorr_long(j,i)+
& wcorr5*gradcorr5_long(j,i)+
& +gradafm(j,i)
& +welec*gshieldc(j,i)
& +welec*gshieldc_loc(j,i)
+ & +wcorr*gshieldc_ec(j,i)
+ & +wcorr*gshieldc_loc_ec(j,i)
+ & +wturn3*gshieldc_t3(j,i)
+ & +wturn3*gshieldc_loc_t3(j,i)
+ & +wturn4*gshieldc_t4(j,i)
+ & +wturn4*gshieldc_loc_t4(j,i)
+ & +wel_loc*gshieldc_ll(j,i)
+ & +wel_loc*gshieldc_loc_ll(j,i)
+
+
+
#endif
& +wscloc*gsclocx(j,i)
& +wliptran*gliptranx(j,i)
& +welec*gshieldx(j,i)
+ & +wcorr*gshieldx_ec(j,i)
+ & +wturn3*gshieldx_t3(j,i)
+ & +wturn4*gshieldx_t4(j,i)
+ & +wel_loc*gshieldx_ll(j,i)
+
+
+
enddo
enddo
#ifdef DEBUG
if (shield_mode.gt.0) then
C fac_shield(i)=0.4
C fac_shield(j)=0.6
- el1=el1*fac_shield(i)*fac_shield(j)
- el2=el2*fac_shield(i)*fac_shield(j)
+ el1=el1*fac_shield(i)**2*fac_shield(j)**2
+ el2=el2*fac_shield(i)**2*fac_shield(j)**2
eesij=(el1+el2)
ees=ees+eesij
else
iresshield=shield_list(ilist,i)
do k=1,3
rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
+ & *2.0
gshieldx(k,iresshield)=gshieldx(k,iresshield)+
& rlocshield
- & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+ & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
iresshield=shield_list(ilist,j)
do k=1,3
rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
+ & *2.0
gshieldx(k,iresshield)=gshieldx(k,iresshield)+
& rlocshield
- & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+ & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
do k=1,3
gshieldc(k,i)=gshieldc(k,i)+
- & grad_shield(k,i)*eesij/fac_shield(i)
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0
gshieldc(k,j)=gshieldc(k,j)+
- & grad_shield(k,j)*eesij/fac_shield(j)
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0
gshieldc(k,i-1)=gshieldc(k,i-1)+
- & grad_shield(k,i)*eesij/fac_shield(i)
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0
gshieldc(k,j-1)=gshieldc(k,j-1)+
- & grad_shield(k,j)*eesij/fac_shield(j)
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0
enddo
endif
cd & (dcosg(k),k=1,3)
do k=1,3
ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
- & fac_shield(i)*fac_shield(j)
+ & fac_shield(i)**2*fac_shield(j)**2
enddo
c do k=1,3
c ghalf=0.5D0*ggg(k)
gelc(k,i)=gelc(k,i)
& +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)**2*fac_shield(j)**2
gelc(k,j)=gelc(k,j)
& +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)**2*fac_shield(j)**2
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
C Contribution to the local-electrostatic energy coming from the i-j pair
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
& +a33*muij(4)
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+C else
+C fac_shield(i)=0.4
+C fac_shield(j)=0.6
+ endif
+ eel_loc_ij=eel_loc_ij
+ & *fac_shield(i)*fac_shield(j)
+C Now derivative over eel_loc
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
+ & /fac_shield(i)
+C & *2.0
+ gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
+ gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
+ & /fac_shield(j)
+C & *2.0
+ gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
+ gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
+ & +rlocshield
+
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc_ll(k,i)=gshieldc_ll(k,i)+
+ & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
+ gshieldc_ll(k,j)=gshieldc_ll(k,j)+
+ & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
+ gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
+ & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
+ gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
+ & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
+ enddo
+ endif
+
+
c write (iout,*) 'i',i,' j',j,itype(i),itype(j),
c & ' eel_loc_ij',eel_loc_ij
C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
C Calculate patrial derivative for theta angle
#ifdef NEWCORR
- geel_loc_ij=a22*gmuij1(1)
+ geel_loc_ij=(a22*gmuij1(1)
& +a23*gmuij1(2)
& +a32*gmuij1(3)
- & +a33*gmuij1(4)
+ & +a33*gmuij1(4))
+ & *fac_shield(i)*fac_shield(j)
c write(iout,*) "derivative over thatai"
c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
c & a33*gmuij1(4)
& +a33*gmuij2(4)
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
& geel_loc_ij*wel_loc
+ & *fac_shield(i)*fac_shield(j)
+
c Derivative over j residue
geel_loc_ji=a22*gmuji1(1)
& +a23*gmuji1(2)
gloc(nphi+j,icg)=gloc(nphi+j,icg)+
& geel_loc_ji*wel_loc
+ & *fac_shield(i)*fac_shield(j)
+
geel_loc_ji=
& +a22*gmuji2(1)
& +a23*gmuji2(2)
c & a33*gmuji2(4)
gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
& geel_loc_ji*wel_loc
+ & *fac_shield(i)*fac_shield(j)
#endif
cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
C Partial derivatives in virtual-bond dihedral angles gamma
if (i.gt.1)
& gel_loc_loc(i-1)=gel_loc_loc(i-1)+
- & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
- & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
+ & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
+ & *fac_shield(i)*fac_shield(j)
+
gel_loc_loc(j-1)=gel_loc_loc(j-1)+
- & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
- & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
+ & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
+ & *fac_shield(i)*fac_shield(j)
C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
do l=1,3
- ggg(l)=agg(l,1)*muij(1)+
- & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
+ ggg(l)=(agg(l,1)*muij(1)+
+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
cgrad ghalf=0.5d0*ggg(l)
do l=1,3
gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
& aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+
gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
& aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+
gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
& aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+
gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
& aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+
enddo
ENDIF
C Change 12/26/95 to calculate four-body contributions to H-bonding energy
ees0mij=0
endif
c ees0mij=0.0D0
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0d0
+ fac_shield(j)=1.0d0
+ else
+ ees0plist(num_conti,i)=j
+C fac_shield(i)=0.4d0
+C fac_shield(j)=0.6d0
+ endif
ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+ & *fac_shield(i)*fac_shield(j)
ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+ & *fac_shield(i)*fac_shield(j)
C Diagnostics. Comment out or remove after debugging!
c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
gacontp_hb1(k,num_conti,i)=!ghalfp
& +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ & *fac_shield(i)*fac_shield(j)
+
gacontp_hb2(k,num_conti,i)=!ghalfp
& +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ & *fac_shield(i)*fac_shield(j)
+
gacontp_hb3(k,num_conti,i)=gggp(k)
+ & *fac_shield(i)*fac_shield(j)
+
gacontm_hb1(k,num_conti,i)=!ghalfm
& +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ & *fac_shield(i)*fac_shield(j)
+
gacontm_hb2(k,num_conti,i)=!ghalfm
& +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ & *fac_shield(i)*fac_shield(j)
+
gacontm_hb3(k,num_conti,i)=gggm(k)
+ & *fac_shield(i)*fac_shield(j)
+
enddo
C Diagnostics. Comment out or remove after debugging!
cdiag do k=1,3
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.CONTROL'
+ include 'COMMON.SHIELD'
dimension ggg(3)
double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
& e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+C else
+C fac_shield(i)=0.4
+C fac_shield(j)=0.6
+ endif
eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+ & *fac_shield(i)*fac_shield(j)
+ eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
+ & *fac_shield(i)*fac_shield(j)
C Derivatives in theta
gloc(nphi+i,icg)=gloc(nphi+i,icg)
& +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
+ & *fac_shield(i)*fac_shield(j)
gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
& +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
+ & *fac_shield(i)*fac_shield(j)
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
+
+C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
+C Derivatives in shield mode
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eello_t3/fac_shield(i)
+C & *2.0
+ gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eello_t3/fac_shield(i)
+ gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eello_t3/fac_shield(j)
+C & *2.0
+ gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eello_t3/fac_shield(j)
+ gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
+ & +rlocshield
+
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc_t3(k,i)=gshieldc_t3(k,i)+
+ & grad_shield(k,i)*eello_t3/fac_shield(i)
+ gshieldc_t3(k,j)=gshieldc_t3(k,j)+
+ & grad_shield(k,j)*eello_t3/fac_shield(j)
+ gshieldc_t3(k,i-1)=gshieldc_t3(k,i-1)+
+ & grad_shield(k,i)*eello_t3/fac_shield(i)
+ gshieldc_t3(k,j-1)=gshieldc_t3(k,j-1)+
+ & grad_shield(k,j)*eello_t3/fac_shield(j)
+ enddo
+ endif
+
+C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
cd write (2,*) 'i,',i,' j',j,'eello_turn3',
cd & 0.5d0*(pizda(1,1)+pizda(2,2)),
cd & ' eello_turn3_num',4*eello_turn3_num
call transpose2(auxmat2(1,1),auxmat3(1,1))
call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))
+ & *fac_shield(i)*fac_shield(j)
C Derivatives in gamma(i+1)
call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1))
call transpose2(auxmat2(1,1),auxmat3(1,1))
call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
gel_loc_turn3(i+1)=gel_loc_turn3(i+1)
& +0.5d0*(pizda(1,1)+pizda(2,2))
+ & *fac_shield(i)*fac_shield(j)
C Cartesian derivatives
do l=1,3
c ghalf1=0.5d0*agg(l,1)
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
gcorr3_turn(l,i)=gcorr3_turn(l,i)
& +0.5d0*(pizda(1,1)+pizda(2,2))
+ & *fac_shield(i)*fac_shield(j)
+
a_temp(1,1)=aggi1(l,1)!+agg(l,1)
a_temp(1,2)=aggi1(l,2)!+agg(l,2)
a_temp(2,1)=aggi1(l,3)!+agg(l,3)
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1)
& +0.5d0*(pizda(1,1)+pizda(2,2))
+ & *fac_shield(i)*fac_shield(j)
a_temp(1,1)=aggj(l,1)!+ghalf1
a_temp(1,2)=aggj(l,2)!+ghalf2
a_temp(2,1)=aggj(l,3)!+ghalf3
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
gcorr3_turn(l,j)=gcorr3_turn(l,j)
& +0.5d0*(pizda(1,1)+pizda(2,2))
+ & *fac_shield(i)*fac_shield(j)
a_temp(1,1)=aggj1(l,1)
a_temp(1,2)=aggj1(l,2)
a_temp(2,1)=aggj1(l,3)
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
gcorr3_turn(l,j1)=gcorr3_turn(l,j1)
& +0.5d0*(pizda(1,1)+pizda(2,2))
+ & *fac_shield(i)*fac_shield(j)
enddo
return
end
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.CONTROL'
+ include 'COMMON.SHIELD'
dimension ggg(3)
double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
& e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2))
gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2))
gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2))
-
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+C else
+C fac_shield(i)=0.6
+C fac_shield(j)=0.4
+ endif
eello_turn4=eello_turn4-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
+ eello_t4=-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2)
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)')
& 'eturn4',i,j,-(s1+s2+s3),s1,s2,s3
+C Now derivative over shield:
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i)
+C & *2.0
+ gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i)
+ gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j)
+C & *2.0
+ gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j)
+ gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
+ & +rlocshield
+
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc_t4(k,i)=gshieldc_t4(k,i)+
+ & grad_shield(k,i)*eello_t4/fac_shield(i)
+ gshieldc_t4(k,j)=gshieldc_t4(k,j)+
+ & grad_shield(k,j)*eello_t4/fac_shield(j)
+ gshieldc_t4(k,i-1)=gshieldc_t4(k,i-1)+
+ & grad_shield(k,i)*eello_t4/fac_shield(i)
+ gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+
+ & grad_shield(k,j)*eello_t4/fac_shield(j)
+ enddo
+ endif
+
+
+
+
+
+
cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
cd & ' eello_turn4_num',8*eello_turn4_num
#ifdef NEWCORR
gloc(nphi+i,icg)=gloc(nphi+i,icg)
& -(gs13+gsE13+gsEE1)*wturn4
+ & *fac_shield(i)*fac_shield(j)
gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)
& -(gs23+gs21+gsEE2)*wturn4
+ & *fac_shield(i)*fac_shield(j)
+
gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)
& -(gs32+gsE31+gsEE3)*wturn4
+ & *fac_shield(i)*fac_shield(j)
+
c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
c & gs2
#endif
call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3)
+ & *fac_shield(i)*fac_shield(j)
C Derivatives in gamma(i+1)
call transpose2(EUgder(1,1,i+2),e2tder(1,1))
call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1))
call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3)
+ & *fac_shield(i)*fac_shield(j)
C Derivatives in gamma(i+2)
call transpose2(EUgder(1,1,i+3),e3tder(1,1))
call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1))
call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
C Cartesian derivatives
C Derivatives of this turn contributions in DC(i+2)
if (j.lt.nres-1) then
s3=0.5d0*(pizda(1,1)+pizda(2,2))
ggg(l)=-(s1+s2+s3)
gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
enddo
endif
C Remaining derivatives of this turn contribution
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
a_temp(1,1)=aggi1(l,1)
a_temp(1,2)=aggi1(l,2)
a_temp(2,1)=aggi1(l,3)
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
a_temp(1,1)=aggj(l,1)
a_temp(1,2)=aggj(l,2)
a_temp(2,1)=aggj(l,3)
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
a_temp(1,1)=aggj1(l,1)
a_temp(1,2)=aggj1(l,2)
a_temp(2,1)=aggj1(l,3)
s3=0.5d0*(pizda(1,1)+pizda(2,2))
c write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3
gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
enddo
return
end
return
end
#endif
+C----------------------------------------------------------------------------------
+C The rigorous attempt to derive energy function
+ subroutine etor_kcc(etors,edihcnstr)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.TORSION'
+ include 'COMMON.INTERACT'
+ include 'COMMON.DERIV'
+ include 'COMMON.CHAIN'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.FFIELD'
+ include 'COMMON.TORCNSTR'
+ include 'COMMON.CONTROL'
+ logical lprn
+ double precision thybt1(maxtermkcc),thybt2(maxtermkcc)
+C Set lprn=.true. for debugging
+ lprn=.false.
+c lprn=.true.
+C print *,"wchodze kcc"
+ if (tor_mode.ne.2) then
+ etors=0.0D0
+ endif
+ do i=iphi_start,iphi_end
+C ANY TWO ARE DUMMY ATOMS in row CYCLE
+c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
+c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or.
+c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
+ if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+ & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+ itori=itortyp_kcc(itype(i-2))
+ itori1=itortyp_kcc(itype(i-1))
+ phii=phi(i)
+ glocig=0.0D0
+ glocit1=0.0d0
+ glocit2=0.0d0
+ sumnonchebyshev=0.0d0
+ sumchebyshev=0.0d0
+C to avoid multiple devision by 2
+ theti22=0.5d0*theta(i)
+C theta 12 is the theta_1 /2
+C theta 22 is theta_2 /2
+ theti12=0.5d0*theta(i-1)
+C and appropriate sinus function
+ sinthet2=dsin(theta(i))
+ sinthet1=dsin(theta(i-1))
+ costhet1=dcos(theta(i-1))
+ costhet2=dcos(theta(i))
+C to speed up lets store its mutliplication
+ sint1t2=sinthet2*sinthet1
+C \sum_{i=1}^n (sin(theta_1) * sin(theta_2))^n * (c_n* cos(n*gamma)
+C +d_n*sin(n*gamma)) *
+C \sum_{i=1}^m (1+a_m*Tb_m(cos(theta_1 /2))+b_m*Tb_m(cos(theta_2 /2)))
+C we have two sum 1) Non-Chebyshev which is with n and gamma
+ do j=1,nterm_kcc(itori,itori1)
+
+ v1ij=v1_kcc(j,itori,itori1)
+ v2ij=v2_kcc(j,itori,itori1)
+C v1ij is c_n and d_n in euation above
+ cosphi=dcos(j*phii)
+ sinphi=dsin(j*phii)
+ sint1t2n=sint1t2**j
+ sumnonchebyshev=
+ & sint1t2n*(v1ij*cosphi+v2ij*sinphi)
+ actval=sint1t2n*(v1ij*cosphi+v2ij*sinphi)
+C etors=etors+sint1t2n*(v1ij*cosphi+v2ij*sinphi)
+C if (energy_dec) etors_ii=etors_ii+
+C & v1ij*cosphi+v2ij*sinphi
+C glocig is the gradient local i site in gamma
+ glocig=j*(v2ij*cosphi-v1ij*sinphi)*sint1t2n
+C now gradient over theta_1
+ glocit1=actval/sinthet1*j*costhet1
+ glocit2=actval/sinthet2*j*costhet2
+
+C now the Czebyshev polinominal sum
+ do k=1,nterm_kcc_Tb(itori,itori1)
+ thybt1(k)=v1_chyb(k,j,itori,itori1)
+ thybt2(k)=v2_chyb(k,j,itori,itori1)
+C thybt1(k)=0.0
+C thybt2(k)=0.0
+ enddo
+ sumth1thyb=tschebyshev
+ & (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theti12)**2)
+ gradthybt1=gradtschebyshev
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1),
+ & dcos(theti12)**2)
+ & *dcos(theti12)*(-dsin(theti12))
+ sumth2thyb=tschebyshev
+ & (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theti22)**2)
+ gradthybt2=gradtschebyshev
+ & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),
+ & dcos(theti22)**2)
+ & *dcos(theti22)*(-dsin(theti22))
+C print *, sumth1thyb, gradthybt1, sumth2thyb, gradthybt2,
+C & gradtschebyshev
+C & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),
+C & dcos(theti22)**2),
+C & dsin(theti22)
+
+C now overal sumation
+ etors=etors+(1.0d0+sumth1thyb+sumth2thyb)*sumnonchebyshev
+C print *,"sumnon", sumnonchebyshev,sumth1thyb+sumth2thyb
+C derivative over gamma
+ gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig
+ & *(1.0d0+sumth1thyb+sumth2thyb)
+C derivative over theta1
+ gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*
+ & (glocit1*(1.0d0+sumth1thyb+sumth2thyb)+
+ & sumnonchebyshev*gradthybt1)
+C now derivative over theta2
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*
+ & (glocit2*(1.0d0+sumth1thyb+sumth2thyb)+
+ & sumnonchebyshev*gradthybt2)
+ enddo
+ enddo
+
+C gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
+! 6/20/98 - dihedral angle constraints
+ if (tor_mode.ne.2) then
+ edihcnstr=0.0d0
+c do i=1,ndih_constr
+ do i=idihconstr_start,idihconstr_end
+ itori=idih_constr(i)
+ phii=phi(itori)
+ difi=pinorm(phii-phi0(i))
+ if (difi.gt.drange(i)) then
+ difi=difi-drange(i)
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+ else
+ difi=0.0
+ endif
+ enddo
+ endif
+ return
+ end
+
+C The rigorous attempt to derive energy function
+ subroutine ebend_kcc(etheta,ethetacnstr)
+
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.TORSION'
+ include 'COMMON.INTERACT'
+ include 'COMMON.DERIV'
+ include 'COMMON.CHAIN'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.FFIELD'
+ include 'COMMON.TORCNSTR'
+ include 'COMMON.CONTROL'
+ logical lprn
+ double precision thybt1(maxtermkcc)
+C Set lprn=.true. for debugging
+ lprn=.false.
+c lprn=.true.
+C print *,"wchodze kcc"
+ if (tormode.ne.2) etheta=0.0D0
+ do i=ithet_start,ithet_end
+c print *,i,itype(i-1),itype(i),itype(i-2)
+ if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+ & .or.itype(i).eq.ntyp1) cycle
+ iti=itortyp_kcc(itype(i-1))
+ sinthet=dsin(theta(i)/2.0d0)
+ costhet=dcos(theta(i)/2.0d0)
+ do j=1,nbend_kcc_Tb(iti)
+ thybt1(j)=v1bend_chyb(j,iti)
+ enddo
+ sumth1thyb=tschebyshev
+ & (1,nbend_kcc_Tb(iti),thybt1(1),costhet)
+ ihelp=nbend_kcc_Tb(iti)-1
+ gradthybt1=gradtschebyshev
+ & (0,ihelp,thybt1(1),costhet)
+ etheta=etheta+sumth1thyb
+C print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0)
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*
+ & gradthybt1*sinthet*(-0.5d0)
+ enddo
+ if (tormode.ne.2) then
+ ethetacnstr=0.0d0
+C print *,ithetaconstr_start,ithetaconstr_end,"TU"
+ do i=ithetaconstr_start,ithetaconstr_end
+ itheta=itheta_constr(i)
+ thetiii=theta(itheta)
+ difi=pinorm(thetiii-theta_constr0(i))
+ if (difi.gt.theta_drange(i)) then
+ difi=difi-theta_drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else
+ difi=0.0
+ endif
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+ & i,itheta,rad2deg*thetiii,
+ & rad2deg*theta_constr0(i), rad2deg*theta_drange(i),
+ & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+ & gloc(itheta+nphi-2,icg)
+ endif
+ enddo
+ endif
+ return
+ end
c------------------------------------------------------------------------------
subroutine eback_sc_corr(esccor)
c 7/21/2007 Correlations between the backbone-local and side-chain-local
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.SHIELD'
double precision gx(3),gx1(3)
logical lprn
lprn=.false.
include 'COMMON.CONTACTS'
include 'COMMON.CHAIN'
include 'COMMON.CONTROL'
+ include 'COMMON.SHIELD'
double precision gx(3),gx1(3)
integer num_cont_hb_old(maxres)
logical lprn,ldone
call calc_eello(i,jp,i+1,jp1,jj,kk)
if (wcorr4.gt.0.0d0)
& ecorr=ecorr+eello4(i,jp,i+1,jp1,jj,kk)
+CC & *fac_shield(i)**2*fac_shield(j)**2
if (energy_dec.and.wcorr4.gt.0.0d0)
1 write (iout,'(a6,4i5,0pf7.3)')
2 'ecorr4',i,j,i+1,j1,eello4(i,jp,i+1,jp1,jj,kk)
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.SHIELD'
+ include 'COMMON.CONTROL'
double precision gx(3),gx1(3)
logical lprn
lprn=.false.
+C print *,"wchodze",fac_shield(i),shield_mode
eij=facont_hb(jj,i)
ekl=facont_hb(kk,k)
ees0pij=ees0p(jj,i)
ees0mkl=ees0m(kk,k)
ekont=eij*ekl
ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+C*
+C & fac_shield(i)**2*fac_shield(j)**2
cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
C Following 4 lines for diagnostics.
cd ees0pkl=0.0D0
cgrad enddo
cgrad enddo
c write (iout,*) "ehbcorr",ekont*ees
+C print *,ekont,ees,i,k
ehbcorr=ekont*ees
+C now gradient over shielding
+C return
+ if (shield_mode.gt.0) then
+ j=ees0plist(jj,i)
+ l=ees0plist(kk,k)
+C print *,i,j,fac_shield(i),fac_shield(j),
+C &fac_shield(k),fac_shield(l)
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (fac_shield(k).gt.0).and.(fac_shield(l).gt.0)) then
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do m=1,3
+ rlocshield=grad_shield_side(m,ilist,i)*ehbcorr/fac_shield(i)
+C & *2.0
+ gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(m,ilist,i)*ehbcorr/fac_shield(i)
+ gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+ &+rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do m=1,3
+ rlocshield=grad_shield_side(m,ilist,j)*ehbcorr/fac_shield(j)
+C & *2.0
+ gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(m,ilist,j)*ehbcorr/fac_shield(j)
+ gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+
+ do ilist=1,ishield_list(k)
+ iresshield=shield_list(ilist,k)
+ do m=1,3
+ rlocshield=grad_shield_side(m,ilist,k)*ehbcorr/fac_shield(k)
+C & *2.0
+ gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(m,ilist,k)*ehbcorr/fac_shield(k)
+ gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(l)
+ iresshield=shield_list(ilist,l)
+ do m=1,3
+ rlocshield=grad_shield_side(m,ilist,l)*ehbcorr/fac_shield(l)
+C & *2.0
+ gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(m,ilist,l)*ehbcorr/fac_shield(l)
+ gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+C print *,gshieldx(m,iresshield)
+ do m=1,3
+ gshieldc_ec(m,i)=gshieldc_ec(m,i)+
+ & grad_shield(m,i)*ehbcorr/fac_shield(i)
+ gshieldc_ec(m,j)=gshieldc_ec(m,j)+
+ & grad_shield(m,j)*ehbcorr/fac_shield(j)
+ gshieldc_ec(m,i-1)=gshieldc_ec(m,i-1)+
+ & grad_shield(m,i)*ehbcorr/fac_shield(i)
+ gshieldc_ec(m,j-1)=gshieldc_ec(m,j-1)+
+ & grad_shield(m,j)*ehbcorr/fac_shield(j)
+
+ gshieldc_ec(m,k)=gshieldc_ec(m,k)+
+ & grad_shield(m,k)*ehbcorr/fac_shield(k)
+ gshieldc_ec(m,l)=gshieldc_ec(m,l)+
+ & grad_shield(m,l)*ehbcorr/fac_shield(l)
+ gshieldc_ec(m,k-1)=gshieldc_ec(m,k-1)+
+ & grad_shield(m,k)*ehbcorr/fac_shield(k)
+ gshieldc_ec(m,l-1)=gshieldc_ec(m,l-1)+
+ & grad_shield(m,l)*ehbcorr/fac_shield(l)
+
+ enddo
+ endif
+ endif
return
end
#ifdef MOMENT
enddo
return
end
+C--------------------------------------------------------------------------
+ double precision function tschebyshev(m,n,x,y)
+ implicit none
+ include "DIMENSIONS"
+ integer i,m,n
+ double precision x(n),y,yy(0:maxvar),aux
+c Tschebyshev polynomial. Note that the first term is omitted
+c m=0: the constant term is included
+c m=1: the constant term is not included
+ yy(0)=1.0d0
+ yy(1)=y
+ do i=2,n
+ yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i)*yy(i)
+ enddo
+ tschebyshev=aux
+ return
+ end
+C--------------------------------------------------------------------------
+ double precision function gradtschebyshev(m,n,x,y)
+ implicit none
+ include "DIMENSIONS"
+ integer i,m,n
+ double precision x(n+1),y,yy(0:maxvar),aux
+c Tschebyshev polynomial. Note that the first term is omitted
+c m=0: the constant term is included
+c m=1: the constant term is not included
+ yy(0)=1.0d0
+ yy(1)=2.0d0*y
+ do i=2,n
+ yy(i)=2*y*yy(i-1)-yy(i-2)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i+1)*yy(i)*(i+1)
+C print *, x(i+1),yy(i),i
+ enddo
+ gradtschebyshev=aux
+ return
+ end