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=00862f24900104f5e82d9398b6736fa125d10a1d;hb=ce23246258bf50c33a76d6b2685ff2a3f76f896d;hp=109af351c026cf52a7ab4975e320166eb0a54e1d;hpb=7ac044004baeac41c3c91c38a1b3c3546287e624;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 109af35..00862f2 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -141,7 +141,7 @@ C Introduction of shielding effect first for each peptide group C the shielding factor is set this factor is describing how each C peptide group is shielded by side-chains C the matrix - shield_fac(i) the i index describe the ith between i and i+1 - write (iout,*) "shield_mode",shield_mode +C write (iout,*) "shield_mode",shield_mode if (shield_mode.gt.0) then call set_shield_fac endif @@ -546,6 +546,7 @@ c enddo & wstrain*ghpbc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) + & +welec*gshieldc(j,i) enddo enddo @@ -564,6 +565,7 @@ c enddo & wstrain*ghpbc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) + & +welec*gshieldc(j,i) enddo enddo @@ -682,6 +684,13 @@ c enddo do i=-1,nct do j=1,3 #ifdef SPLITELE +C print *,gradbufc(1,13) +C print *,welec*gelc(1,13) +C print *,wel_loc*gel_loc(1,13) +C print *,0.5d0*(wscp*gvdwc_scpp(1,13)) +C print *,welec*gelc_long(1,13)+wvdwpp*gvdwpp(1,13) +C print *,wel_loc*gel_loc_long(1,13) +C print *,gradafm(1,13),"AFM" gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ & 0.5d0*(wscp*gvdwc_scpp(j,i)+ @@ -702,6 +711,10 @@ c enddo & +wscloc*gscloc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + + #else gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ @@ -723,6 +736,9 @@ c enddo & +wscloc*gscloc(j,i) & +wliptran*gliptranc(j,i) & +gradafm(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + #endif gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ @@ -731,6 +747,7 @@ c enddo & wsccor*gsccorx(j,i) & +wscloc*gsclocx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) enddo enddo #ifdef DEBUG @@ -2763,6 +2780,17 @@ c write(iout,*) 'b1=',b1(1,i-2) c write (iout,*) 'theta=', theta(i-1) enddo #else + if (i.gt. nnt+2 .and. i.lt.nct+2) then + iti = itortyp(itype(i-2)) + else + iti=ntortyp+1 + endif +c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itortyp(itype(i-1)) + else + iti1=ntortyp+1 + endif b1(1,i-2)=b(3,iti) b1(2,i-2)=b(5,iti) b2(1,i-2)=b(2,iti) @@ -2917,6 +2945,7 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo +C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2) c write (iout,*) 'mu ',mu(:,i-2),i-2 cd write (iout,*) 'mu1',mu1(:,i-2) cd write (iout,*) 'mu2',mu2(:,i-2) @@ -3429,7 +3458,9 @@ C do zshift=-1,1 c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c +CTU KURWA do i=iatel_s,iatel_e +C do i=75,75 if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds @@ -3486,7 +3517,9 @@ c endif c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) +C I TU KURWA do j=ielstart(i),ielend(i) +C do j=16,17 C write (iout,*) i,j if (j.le.1) cycle if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 @@ -3669,12 +3702,20 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac C MARYSIA - eesij=(el1+el2) +C eesij=(el1+el2) 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 - ees=ees+eesij*fac_shield(i)*fac_shield(j) +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) + 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 @@ -3700,22 +3741,99 @@ C erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij + * * Radial derivatives. First process both termini of the fragment (i,j) * 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) + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) + 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) + gshieldx(k,iresshield)=gshieldx(k,iresshield)+ + & rlocshield + & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) + gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield + +C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) +C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C if (iresshield.gt.j) then +C do ishi=j+1,iresshield-1 +C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield +C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C +C enddo +C else +C do ishi=iresshield,j +C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield +C & -grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) +C enddo +C endif + enddo + enddo + + do k=1,3 + gshieldc(k,i)=gshieldc(k,i)+ + & grad_shield(k,i)*eesij/fac_shield(i) + gshieldc(k,j)=gshieldc(k,j)+ + & grad_shield(k,j)*eesij/fac_shield(j) + gshieldc(k,i-1)=gshieldc(k,i-1)+ + & grad_shield(k,i)*eesij/fac_shield(i) + gshieldc(k,j-1)=gshieldc(k,j-1)+ + & grad_shield(k,j)*eesij/fac_shield(j) + + enddo + endif c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf c gelc(k,j)=gelc(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end +C print *,"before", gelc_long(1,i), gelc_long(1,j) do k=1,3 gelc_long(k,j)=gelc_long(k,j)+ggg(k) +C & +grad_shield(k,j)*eesij/fac_shield(j) gelc_long(k,i)=gelc_long(k,i)-ggg(k) +C & +grad_shield(k,i)*eesij/fac_shield(i) +C gelc_long(k,i-1)=gelc_long(k,i-1) +C & +grad_shield(k,i)*eesij/fac_shield(i) +C gelc_long(k,j-1)=gelc_long(k,j-1) +C & +grad_shield(k,j)*eesij/fac_shield(j) enddo +C print *,"bafter", gelc_long(1,i), gelc_long(1,j) + * * Loop over residues i+1 thru j-1. * @@ -3764,8 +3882,11 @@ C MARYSIA * Radial derivatives. First process both termini of the fragment (i,j) * ggg(1)=fac*xj +C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j) ggg(2)=fac*yj +C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j) ggg(3)=fac*zj +C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j) c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf @@ -3808,7 +3929,8 @@ c 9/28/08 AL Gradient compotents will be summed only at the end 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) + ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* + & fac_shield(i)*fac_shield(j) enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -3824,16 +3946,21 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo +C print *,"before22", gelc_long(1,i), gelc_long(1,j) do k=1,3 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) + & +((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) 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) + & +((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) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo +C print *,"before33", gelc_long(1,i), gelc_long(1,j) + C MARYSIA c endif !sscale IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 @@ -4036,7 +4163,7 @@ C Contribution to the local-electrostatic energy coming from the i-j pair & +a33*muij(4) c write (iout,*) 'i',i,' j',j,itype(i),itype(j), c & ' eel_loc_ij',eel_loc_ij -c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4) +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) @@ -5783,6 +5910,7 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end + if (i.le.2) cycle 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 @@ -5799,6 +5927,15 @@ C print *,i,theta(i) sinkt(k)=dsin(k*theti2) enddo C print *,ethetai + if (i.eq.3) then + phii=0.0d0 + ityp1=nthetyp+1 + do k=1,nsingle + cosph1(k)=0.0d0 + sinph1(k)=0.0d0 + enddo + else + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) @@ -5820,6 +5957,7 @@ C propagation of chirality for glycine type sinph1(k)=0.0d0 enddo endif + endif if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) @@ -10530,11 +10668,12 @@ C first for shielding is setting of function of side-chains 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/ + &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) + &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 @@ -10543,14 +10682,17 @@ 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(k+nres,j)-(c(i,j)+c(i+1,j))/2.0d0 + 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(k+nres,j)-c(k,j) - pept_group(j)=c(i,j)-c(i+1,j) + 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 @@ -10558,8 +10700,14 @@ C lets have their lenght 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 @@ -10567,7 +10715,7 @@ 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))=k + 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 @@ -10577,20 +10725,29 @@ C Lets have the sscale value else scale_fac_dist=-sh_frac_dist*sh_frac_dist & *(2.0*sh_frac_dist-3.0d0) - fac_help_scale=6.0*(scale_fac_dist-scale_fac_dist**2) + 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+short**2/dist_pep_side**2) + costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2) C now costhet_grad - costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**3 +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 @@ -10602,11 +10759,30 @@ C pep_side0pept_group is vector multiplication do j=1,3 pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) enddo - fac_alfa_sin=1.0-(pep_side0pept_group/ - & (dist_pep_side*dist_side_calf))**2 + 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 - cosphi=1.0d0/dsqrt(1+rkprim**2/dist_pep_side**2) +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 C now the gradient... @@ -10614,27 +10790,29 @@ C grad_shield is gradient of Calfa for peptide groups do j=1,3 grad_shield(j,i)=grad_shield(j,i) C gradient po skalowaniu - & +sh_frac_dist_grad(j)*VofOverlap + & +(sh_frac_dist_grad(j) C gradient po costhet - &+scale_fac_dist*costhet_grad(j) + &-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)*VofOverlap*2.0d0 - & +scale_fac_dist*costhet_grad(j)*2.0d0 -C grad_shield_side_ca is Calfa sidechain gradient - grad_shield_side_ca(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 -C if ((cosphi.le.0.0).or.(costhet.le.0.0)) write(iout,*) "ERROR", -C & cosphi,costhet -C now should be fac_side_grad(k) which will be gradient of factor k which also -C affect the gradient of peptide group i fac_pept_grad(i) and i+1 - write(2,*) "myvolume",VofOverlap,VSolvSphere_div,VolumeTotal - enddo -C write(2,*) "TOTAL VOLUME",i,VolumeTotal -C the scaling factor of the shielding effect + enddo fac_shield(i)=VolumeTotal*div77_81+div4_81 - write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) +C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) enddo return end