& +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
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
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)
ees0mij=ees0m(jj,i)
ees0mkl=ees0m(kk,k)
ekont=eij*ekl
- ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)*
- & fac_shield(i)**2*fac_shield(j)**2
+ 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