include 'COMMON.INTERACT'
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
+ include 'COMMON.SHIELD'
+ include 'COMMON.CONTROL'
double precision fact(6)
cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot
cd print *,'nnt=',nnt,' nct=',nct
C
C Calculate electrostatic (H-bonding) energy of the main chain.
C
- 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+ 106 continue
+ if (shield_mode.eq.1) then
+ call set_shield_fac
+ else if (shield_mode.eq.2) then
+ call set_shield_fac2
+ endif
+ call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
C write(iout,*) 'po eelec'
C Calculate excluded-volume interaction energy between peptide groups
call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1)
c write (*,*) 'n_corr=',n_corr,' n_corr1=',n_corr1
c print *,ecorr,ecorr5,ecorr6,eturn6
+ else
+ ecorr=0.0d0
+ ecorr5=0.0d0
+ ecorr6=0.0d0
+ eturn6=0.0d0
endif
if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
endif
c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
#ifdef SPLITELE
+ if (shield_mode.gt.0) then
+ etot=fact(1)*wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2
+ & +welec*fact(1)*ees
+ & +fact(1)*wvdwpp*evdw1
+ & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
+ & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
+ & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
+ & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+ & +wliptran*eliptran
+ else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
& +wvdwpp*evdw1
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
& +wliptran*eliptran
+ endif
#else
+ if (shield_mode.gt.0) then
+ etot=fact(1)wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2
+ & +welec*fact(1)*(ees+evdw1)
+ & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
+ & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
+ & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
+ & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+ & +wliptran*eliptran
+ else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
& +wliptran*eliptran
+ endif
#endif
energia(0)=etot
energia(1)=evdw
#ifdef SPLITELE
do i=1,nct
do j=1,3
+ if (shield_mode.eq.0) then
gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+
& welec*fact(1)*gelc(j,i)+wvdwpp*gvdwpp(j,i)+
& wbond*gradb(j,i)+
& wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
& wsccor*fact(2)*gsccorx(j,i)
& +wliptran*gliptranx(j,i)
+ else
+ gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)
+ & +fact(1)*wscp*gvdwc_scp(j,i)+
+ & welec*fact(1)*gelc(j,i)+fact(1)*wvdwpp*gvdwpp(j,i)+
+ & wbond*gradb(j,i)+
+ & wstrain*ghpbc(j,i)+
+ & wcorr*fact(3)*gradcorr(j,i)+
+ & wel_loc*fact(2)*gel_loc(j,i)+
+ & wturn3*fact(2)*gcorr3_turn(j,i)+
+ & wturn4*fact(3)*gcorr4_turn(j,i)+
+ & wcorr5*fact(4)*gradcorr5(j,i)+
+ & wcorr6*fact(5)*gradcorr6(j,i)+
+ & wturn6*fact(5)*gcorr6_turn(j,i)+
+ & wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(j,i)
+ gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i)
+ & +fact(1)*wscp*gradx_scp(j,i)+
+ & wbond*gradbx(j,i)+
+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
+ & wsccor*fact(2)*gsccorx(j,i)
+ & +wliptran*gliptranx(j,i)
+
+ endif
enddo
#else
do i=1,nct
do j=1,3
+ if (shield_mode.eq.0) then
gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+
& welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+
& wbond*gradb(j,i)+
& wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
& wsccor*fact(1)*gsccorx(j,i)
& +wliptran*gliptranx(j,i)
+ else
+ gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)+
+ & fact(1)*wscp*gvdwc_scp(j,i)+
+ & welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+
+ & wbond*gradb(j,i)+
+ & wcorr*fact(3)*gradcorr(j,i)+
+ & wel_loc*fact(2)*gel_loc(j,i)+
+ & wturn3*fact(2)*gcorr3_turn(j,i)+
+ & wturn4*fact(3)*gcorr4_turn(j,i)+
+ & wcorr5*fact(4)*gradcorr5(j,i)+
+ & wcorr6*fact(5)*gradcorr6(j,i)+
+ & wturn6*fact(5)*gcorr6_turn(j,i)+
+ & wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(j,i)
+ gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i)+
+ & fact(1)*wscp*gradx_scp(j,i)+
+ & wbond*gradbx(j,i)+
+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
+ & wsccor*fact(1)*gsccorx(j,i)
+ & +wliptran*gliptranx(j,i)
+ endif
enddo
#endif
enddo
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
+ include 'COMMON.SHIELD'
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
c write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij
C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
+ if (shield_mode.gt.0) then
+C fac_shield(i)=0.4
+C fac_shield(j)=0.6
+ 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
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ eesij=(el1+el2)
+ ees=ees+eesij
+ endif
evdw1=evdw1+evdwij*sss
c write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)')
c &'evdw1',i,j,evdwij
ggg(1)=facel*xj
ggg(2)=facel*yj
ggg(3)=facel*zj
+ 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)*eesij/fac_shield(i)
+ & *2.0
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +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)
+C if (iresshield.gt.i) then
+C do ishi=i+1,iresshield-1
+C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C else
+C do ishi=iresshield,i
+C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C endif
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ 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)*2.0
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc(k,i)=gshieldc(k,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)*2.0
+ gshieldc(k,i-1)=gshieldc(k,i-1)+
+ & 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)*2.0
+
+ enddo
+ endif
+
do k=1,3
ghalf=0.5D0*ggg(k)
gelc(k,i)=gelc(k,i)+ghalf
cd & (dcosg(k),k=1,3)
do k=1,3
ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ & *fac_shield(i)**2*fac_shield(j)**2
enddo
do k=1,3
ghalf=0.5D0*ggg(k)
gelc(k,i)=gelc(k,i)+ghalf
& +(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)**2*fac_shield(j)**2
+
gelc(k,j)=gelc(k,j)+ghalf
& +(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)**2*fac_shield(j)**2
enddo
do k=i+1,j-1
do l=1,3
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 write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
C write (iout,'(a6,2i5,0pf7.3)')
C & 'eelloc',i,j,eel_loc_ij
C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
+C eel_loc=eel_loc+eel_loc_ij
+ 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
eel_loc=eel_loc+eel_loc_ij
+
C Partial derivatives in virtual-bond dihedral angles gamma
if (calc_grad) then
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)
+
cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
cd write(iout,*) 'agg ',agg
cd write(iout,*) 'aggi ',aggi
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)
+
enddo
do k=i+2,j2
do l=1,3
enddo
C Remaining derivatives of eello
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)
- 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)
- 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)
- 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)
+ 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
ENDIF
ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
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)
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
endif
C Diagnostics. Comment out or remove after debugging!
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
+ include 'COMMON.SHIELD'
+ include 'COMMON.CONTROL'
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(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
call transpose2(auxmat(1,1),auxmat1(1,1))
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(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)
+
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
if (calc_grad) then
+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 Derivatives in gamma(i)
call matmat2(EUgder(1,1,i+1),EUg(1,1,i+2),auxmat2(1,1))
call transpose2(auxmat2(1,1),pizda(1,1))
call matmat2(a_temp(1,1),pizda(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
a_temp(1,1)=aggi(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)
a_temp(1,2)=aggi1(l,2)
a_temp(2,1)=aggi1(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)
a_temp(1,2)=aggj(l,2)
a_temp(2,1)=aggj(l,3)
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
endif
else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(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.4
+C fac_shield(j)=0.6
+ 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)
+
cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
cd & ' eello_turn4_num',8*eello_turn4_num
C Derivatives in gamma(i)
if (calc_grad) then
+ 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
call transpose2(EUgder(1,1,i+1),e1tder(1,1))
call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
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(auxmat(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
do l=1,3
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)
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
+
enddo
endif
178 continue
C Set lprn=.true. for debugging
lprn=.false.
eturn6=0.0d0
+ ecorr6=0.0d0
#ifdef MPL
n_corr=0
n_corr1=0
cd write (iout,*) '******eturn6: i,j,i+1,j1',i,j,i+1,j1
eturn6=eturn6+eello_turn6(i,jj,kk)
cd write (2,*) 'multibody_eello:eturn6',eturn6
+ else if ((wturn6.eq.0.0d0).and.(wcorr6.eq.0.0d0)) then
+ eturn6=0.0d0
+ ecorr6=0.0d0
endif
+
ENDIF
1111 continue
else if (j1.eq.j) then
enddo ! kk
enddo ! jj
enddo ! i
+ write (iout,*) "eturn6",eturn6,ecorr6
return
end
c------------------------------------------------------------------------------
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SHIELD'
double precision gx(3),gx1(3)
logical lprn
lprn=.false.
c & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
C Calculate the multi-body contribution to energy.
- ecorr=ecorr+ekont*ees
+C ecorr=ecorr+ekont*ees
if (calc_grad) then
C Calculate multi-body contributions to the gradient.
do ll=1,3
& ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+
& coeffm*ees0mij*gacontm_hb3(ll,kk,k))
enddo
- enddo
+ enddo
+ 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
endif
ehbcorr=ekont*ees
return
end
C-----------------------------------------------------------------------
+ subroutine set_shield_fac
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SHIELD'
+ include 'COMMON.INTERACT'
+C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
+ double precision div77_81/0.974996043d0/,
+ &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
+
+C the vector between center of side_chain and peptide group
+ double precision pep_side(3),long,side_calf(3),
+ &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+ &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C the line belowe needs to be changed for FGPROC>1
+ do i=1,nres-1
+ if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ ishield_list(i)=0
+Cif there two consequtive dummy atoms there is no peptide group between them
+C the line below has to be changed for FGPROC>1
+ VolumeTotal=0.0
+ do k=1,nres
+ if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+ dist_pep_side=0.0
+ dist_side_calf=0.0
+ do j=1,3
+C first lets set vector conecting the ithe side-chain with kth side-chain
+ pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C pep_side(j)=2.0d0
+C and vector conecting the side-chain with its proper calfa
+ side_calf(j)=c(j,k+nres)-c(j,k)
+C side_calf(j)=2.0d0
+ pept_group(j)=c(j,i)-c(j,i+1)
+C lets have their lenght
+ dist_pep_side=pep_side(j)**2+dist_pep_side
+ dist_side_calf=dist_side_calf+side_calf(j)**2
+ dist_pept_group=dist_pept_group+pept_group(j)**2
+ enddo
+ dist_pep_side=dsqrt(dist_pep_side)
+ dist_pept_group=dsqrt(dist_pept_group)
+ dist_side_calf=dsqrt(dist_side_calf)
+ do j=1,3
+ pep_side_norm(j)=pep_side(j)/dist_pep_side
+ side_calf_norm(j)=dist_side_calf
+ enddo
+C now sscale fraction
+ sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C print *,buff_shield,"buff"
+C now sscale
+ if (sh_frac_dist.le.0.0) cycle
+C If we reach here it means that this side chain reaches the shielding sphere
+C Lets add him to the list for gradient
+ ishield_list(i)=ishield_list(i)+1
+C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+C this list is essential otherwise problem would be O3
+ shield_list(ishield_list(i),i)=k
+C Lets have the sscale value
+ if (sh_frac_dist.gt.1.0) then
+ scale_fac_dist=1.0d0
+ do j=1,3
+ sh_frac_dist_grad(j)=0.0d0
+ enddo
+ else
+ scale_fac_dist=-sh_frac_dist*sh_frac_dist
+ & *(2.0*sh_frac_dist-3.0d0)
+ fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2)
+ & /dist_pep_side/buff_shield*0.5
+C remember for the final gradient multiply sh_frac_dist_grad(j)
+C for side_chain by factor -2 !
+ do j=1,3
+ sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C print *,"jestem",scale_fac_dist,fac_help_scale,
+C & sh_frac_dist_grad(j)
+ enddo
+ endif
+C if ((i.eq.3).and.(k.eq.2)) then
+C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist
+C & ,"TU"
+C endif
+
+C this is what is now we have the distance scaling now volume...
+ short=short_r_sidechain(itype(k))
+ long=long_r_sidechain(itype(k))
+ costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2)
+C now costhet_grad
+C costhet=0.0d0
+ costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4
+C costhet_fac=0.0d0
+ do j=1,3
+ costhet_grad(j)=costhet_fac*pep_side(j)
+ enddo
+C remember for the final gradient multiply costhet_grad(j)
+C for side_chain by factor -2 !
+C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+C pep_side0pept_group is vector multiplication
+ pep_side0pept_group=0.0
+ do j=1,3
+ pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+ enddo
+ cosalfa=(pep_side0pept_group/
+ & (dist_pep_side*dist_side_calf))
+ fac_alfa_sin=1.0-cosalfa**2
+ fac_alfa_sin=dsqrt(fac_alfa_sin)
+ rkprim=fac_alfa_sin*(long-short)+short
+C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2)
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4
+
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+ &+cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa/
+ &((dist_pep_side*dist_side_calf))*
+ &((side_calf(j))-cosalfa*
+ &((pep_side(j)/dist_pep_side)*dist_side_calf))
+
+ cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa
+ &/((dist_pep_side*dist_side_calf))*
+ &(pep_side(j)-
+ &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+ enddo
+
+ VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
+ & /VSolvSphere_div
+ & *wshield
+C now the gradient...
+C grad_shield is gradient of Calfa for peptide groups
+C write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist,
+C & costhet,cosphi
+C write(iout,*) "cosphi_compon",i,k,pep_side0pept_group,
+C & dist_pep_side,dist_side_calf,c(1,k+nres),c(1,k),itype(k)
+ do j=1,3
+ grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+ & +(sh_frac_dist_grad(j)
+C gradient po costhet
+ &-scale_fac_dist*costhet_grad(j)/(1.0-costhet)
+ &-scale_fac_dist*(cosphi_grad_long(j))
+ &/(1.0-cosphi) )*div77_81
+ &*VofOverlap
+C grad_shield_side is Cbeta sidechain gradient
+ grad_shield_side(j,ishield_list(i),i)=
+ & (sh_frac_dist_grad(j)*-2.0d0
+ & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
+ & +scale_fac_dist*(cosphi_grad_long(j))
+ & *2.0d0/(1.0-cosphi))
+ & *div77_81*VofOverlap
+
+ grad_shield_loc(j,ishield_list(i),i)=
+ & scale_fac_dist*cosphi_grad_loc(j)
+ & *2.0d0/(1.0-cosphi)
+ & *div77_81*VofOverlap
+ enddo
+ VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+ enddo
+ fac_shield(i)=VolumeTotal*div77_81+div4_81
+C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+ enddo
+ return
+ end
+C--------------------------------------------------------------------------
+C first for shielding is setting of function of side-chains
+ subroutine set_shield_fac2
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SHIELD'
+ include 'COMMON.INTERACT'
+C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
+ double precision div77_81/0.974996043d0/,
+ &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
+
+C the vector between center of side_chain and peptide group
+ double precision pep_side(3),long,side_calf(3),
+ &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+ &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C the line belowe needs to be changed for FGPROC>1
+ do i=1,nres-1
+ if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ ishield_list(i)=0
+Cif there two consequtive dummy atoms there is no peptide group between them
+C the line below has to be changed for FGPROC>1
+ VolumeTotal=0.0
+ do k=1,nres
+ if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+ dist_pep_side=0.0
+ dist_side_calf=0.0
+ do j=1,3
+C first lets set vector conecting the ithe side-chain with kth side-chain
+ pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C pep_side(j)=2.0d0
+C and vector conecting the side-chain with its proper calfa
+ side_calf(j)=c(j,k+nres)-c(j,k)
+C side_calf(j)=2.0d0
+ pept_group(j)=c(j,i)-c(j,i+1)
+C lets have their lenght
+ dist_pep_side=pep_side(j)**2+dist_pep_side
+ dist_side_calf=dist_side_calf+side_calf(j)**2
+ dist_pept_group=dist_pept_group+pept_group(j)**2
+ enddo
+ dist_pep_side=dsqrt(dist_pep_side)
+ dist_pept_group=dsqrt(dist_pept_group)
+ dist_side_calf=dsqrt(dist_side_calf)
+ do j=1,3
+ pep_side_norm(j)=pep_side(j)/dist_pep_side
+ side_calf_norm(j)=dist_side_calf
+ enddo
+C now sscale fraction
+ sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C print *,buff_shield,"buff"
+C now sscale
+ if (sh_frac_dist.le.0.0) cycle
+C If we reach here it means that this side chain reaches the shielding sphere
+C Lets add him to the list for gradient
+ ishield_list(i)=ishield_list(i)+1
+C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+C this list is essential otherwise problem would be O3
+ shield_list(ishield_list(i),i)=k
+C Lets have the sscale value
+ if (sh_frac_dist.gt.1.0) then
+ scale_fac_dist=1.0d0
+ do j=1,3
+ sh_frac_dist_grad(j)=0.0d0
+ enddo
+ else
+ scale_fac_dist=-sh_frac_dist*sh_frac_dist
+ & *(2.0d0*sh_frac_dist-3.0d0)
+ fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2)
+ & /dist_pep_side/buff_shield*0.5d0
+C remember for the final gradient multiply sh_frac_dist_grad(j)
+C for side_chain by factor -2 !
+ do j=1,3
+ sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C sh_frac_dist_grad(j)=0.0d0
+C scale_fac_dist=1.0d0
+C print *,"jestem",scale_fac_dist,fac_help_scale,
+C & sh_frac_dist_grad(j)
+ enddo
+ endif
+C this is what is now we have the distance scaling now volume...
+ short=short_r_sidechain(itype(k))
+ long=long_r_sidechain(itype(k))
+ costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
+ sinthet=short/dist_pep_side*costhet
+C now costhet_grad
+C costhet=0.6d0
+C sinthet=0.8
+ costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4
+C sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet
+C & -short/dist_pep_side**2/costhet)
+C costhet_fac=0.0d0
+ do j=1,3
+ costhet_grad(j)=costhet_fac*pep_side(j)
+ enddo
+C remember for the final gradient multiply costhet_grad(j)
+C for side_chain by factor -2 !
+C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+C pep_side0pept_group is vector multiplication
+ pep_side0pept_group=0.0d0
+ do j=1,3
+ pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+ enddo
+ cosalfa=(pep_side0pept_group/
+ & (dist_pep_side*dist_side_calf))
+ fac_alfa_sin=1.0d0-cosalfa**2
+ fac_alfa_sin=dsqrt(fac_alfa_sin)
+ rkprim=fac_alfa_sin*(long-short)+short
+C rkprim=short
+
+C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2)
+C cosphi=0.6
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4
+ sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/
+ & dist_pep_side**2)
+C sinphi=0.8
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+ &+cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa/
+ &((dist_pep_side*dist_side_calf))*
+ &((side_calf(j))-cosalfa*
+ &((pep_side(j)/dist_pep_side)*dist_side_calf))
+C cosphi_grad_long(j)=0.0d0
+ cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa
+ &/((dist_pep_side*dist_side_calf))*
+ &(pep_side(j)-
+ &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+C cosphi_grad_loc(j)=0.0d0
+ enddo
+C print *,sinphi,sinthet
+ VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet))
+ & /VSolvSphere_div
+C & *wshield
+C now the gradient...
+ do j=1,3
+ grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+ & +(sh_frac_dist_grad(j)*VofOverlap
+C gradient po costhet
+ & +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0*
+ &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+ & sinphi/sinthet*costhet*costhet_grad(j)
+ & +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+ & )*wshield
+C grad_shield_side is Cbeta sidechain gradient
+ grad_shield_side(j,ishield_list(i),i)=
+ & (sh_frac_dist_grad(j)*-2.0d0
+ & *VofOverlap
+ & -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+ &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+ & sinphi/sinthet*costhet*costhet_grad(j)
+ & +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+ & )*wshield
+
+ grad_shield_loc(j,ishield_list(i),i)=
+ & scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+ &(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(
+ & sinthet/sinphi*cosphi*cosphi_grad_loc(j)
+ & ))
+ & *wshield
+ enddo
+ VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+ enddo
+ fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
+C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+C write(2,*) "TU",rpp(1,1),short,long,buff_shield
+ enddo
+ return
+ end
projects / unres.git / blobdiff