X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=a889c75dfd9a6bdde55732701d63794cfe729a6c;hb=9111b59335b350bc16fcac69159ec3d3cbb62ba4;hp=ee55c93e4847c68a7f6e16324c9d41b1aff922b5;hpb=25618f9f83673a7063414fe1e17415d138f58da8;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index ee55c93..a889c75 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -547,6 +547,11 @@ c enddo & +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 @@ -566,6 +571,10 @@ c 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 @@ -713,13 +722,25 @@ C print *,gradafm(1,13),"AFM" & +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)+ @@ -738,6 +759,17 @@ C print *,gradafm(1,13),"AFM" & +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 @@ -748,6 +780,13 @@ C print *,gradafm(1,13),"AFM" & +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 @@ -3708,8 +3747,8 @@ C 12/26/95 - for the evaluation of multi-body H-bonding interactions 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 @@ -3728,7 +3767,8 @@ cd & xmedi,ymedi,zmedi,xj,yj,zj write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') &'evdw1',i,j,evdwij &,iteli,itelj,aaa,evdw1 - write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij + write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, + &fac_shield(i),fac_shield(j) endif C @@ -3755,9 +3795,10 @@ C print *,i,j 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) @@ -3780,9 +3821,10 @@ C endif 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) @@ -3805,13 +3847,13 @@ C endif 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 @@ -3930,7 +3972,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), 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) @@ -3951,11 +3993,11 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) 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 @@ -4161,15 +4203,71 @@ cgrad endif 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) @@ -4185,6 +4283,8 @@ c & a33*gmuij2(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) @@ -4196,6 +4296,8 @@ c & a33*gmuji1(4) 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) @@ -4206,6 +4308,7 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), 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 @@ -4219,15 +4322,19 @@ c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) 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) @@ -4243,12 +4350,20 @@ 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)) + & *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 @@ -4327,8 +4442,18 @@ c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) 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 @@ -4396,17 +4521,29 @@ cgrad ghalfm=0.5D0*gggm(k) 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 @@ -4458,6 +4595,7 @@ C Third- and fourth-order contributions from turns 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), @@ -4498,15 +4636,71 @@ c auxalary matrix for i+2 and constant i+1 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 @@ -4515,12 +4709,14 @@ C Derivatives in gamma(i) 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) @@ -4534,6 +4730,8 @@ c ghalf4=0.5d0*agg(l,4) 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) @@ -4541,6 +4739,7 @@ c ghalf4=0.5d0*agg(l,4) 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 @@ -4548,6 +4747,7 @@ c ghalf4=0.5d0*agg(l,4) 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) @@ -4555,6 +4755,7 @@ c ghalf4=0.5d0*agg(l,4) 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 @@ -4575,6 +4776,7 @@ C Third- and fourth-order contributions from turns 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), @@ -4675,20 +4877,82 @@ c i+3 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 @@ -4704,6 +4968,7 @@ C Derivatives in gamma(i) 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)) @@ -4712,6 +4977,7 @@ C Derivatives in gamma(i+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)) @@ -4723,6 +4989,7 @@ C Derivatives in gamma(i+2) 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 @@ -4742,6 +5009,7 @@ C Derivatives of this turn contributions in DC(i+2) 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 @@ -4760,6 +5028,7 @@ 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) @@ -4774,6 +5043,7 @@ 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+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) @@ -4788,6 +5058,7 @@ 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,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) @@ -4803,6 +5074,7 @@ C Remaining derivatives of this turn contribution 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 @@ -7301,6 +7573,7 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.SHIELD' double precision gx(3),gx1(3) logical lprn lprn=.false. @@ -7719,6 +7992,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding 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 @@ -8021,6 +8295,7 @@ cd write (iout,*) "grij_hb_cont i1",grij_hb_cont(:,jj,i1) 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) @@ -8140,9 +8415,12 @@ c------------------------------------------------------------------------------ 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) @@ -8151,6 +8429,8 @@ c------------------------------------------------------------------------------ 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 @@ -8214,7 +8494,89 @@ cgrad & coeffm*ees0mij*gacontm_hb3(ll,kk,k)) 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 @@ -10776,6 +11138,10 @@ C now costhet_grad & /VSolvSphere_div 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