From: Adam Sieradzan Date: Mon, 30 Apr 2018 07:49:19 +0000 (+0200) Subject: changes in shielding X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;ds=sidebyside;h=7a01d1dcffa9e3d055581b6b4fc937dd118d1841;hp=-c;p=unres4.git changes in shielding --- 7a01d1dcffa9e3d055581b6b4fc937dd118d1841 diff --git a/PARAM/scparm.adam_30_10_2013 b/PARAM/scparm.adam_30_10_2013 index 82dacd1..d5e645c 100644 --- a/PARAM/scparm.adam_30_10_2013 +++ b/PARAM/scparm.adam_30_10_2013 @@ -350,7 +350,7 @@ 0.0347 4.9555 -0.1115 0.4550 -0.8819 0.7342 0.0000 0.0000 0.0000 0.0000 3.9312 7.0874 -0.1922 0.1141 1 0.0000 0.0000 0.0000 0.0000 1.6344 0.0000 0.0000 0.0000 0.5564 0.0531 0.0000 5.0000 0.2736 1.5702 0.0000 0.0000 0.0000 1.1690 0.0000 7.6197 0.1701 1.0347 13.6260 0.0000 1.0000 1.0000 0.0347 4.9555 -0.1115 0.4550 -0.8819 0.7342 0.0000 0.0000 0.0000 0.0000 3.9312 7.0874 -0.1922 0.1141 1 0.0000 0.0000 0.0000 0.0000 1.6344 0.0000 0.0000 0.0000 0.5564 0.0531 0.0000 5.0000 0.2736 1.5702 0.0000 0.0000 0.0000 1.1690 0.0000 7.6197 0.1701 1.0347 13.6260 0.0000 1.0000 1.0000 0.4409 7.0124 -0.7792 -0.9026 -0.9824 0.9981 0.0000 0.0000 0.0000 0.0000 4.2163 6.7741 -0.9979 -0.6323 1 0.0000 0.0000 0.0000 0.0000 1.2736 0.0000 0.0000 0.0000 6.9619 -0.0105 0.0000 5.0000 0.3624 1.4965 0.0000 0.0000 0.0000 1.5081 0.0000 6.4201 0.1437 0.1181 13.7010 0.0000 1.0000 1.0000 - 1.2808 3.6212 0.2182 0.8105 0.9442 -0.8842 0.0000 0.0000 0.0000 0.0000 1.8965 3.4196 0.1846 0.7056 1 0.0000 0.0000 0.0000 0.0000 1.2281 0.0000 0.0000 0.0000 -0.0631 -0.2997 0.0000 5.0000 1.1189 1.3499 0.0000 0.0000 0.0000 1.3003 0.0000 13.2550 11.0080 0.0349 13.4380 0.0000 1.0000 1.0000 !OThr 1/8 epsilon i +1A sigma + 1.2808 2.6212 0.2182 0.8105 0.9442 -0.8842 0.0000 0.0000 0.0000 0.0000 1.8965 3.4196 0.1846 0.7056 1 0.0000 0.0000 0.0000 0.0000 1.2281 0.0000 0.0000 0.0000 -0.0631 -0.2997 0.0000 5.0000 1.1189 1.3499 0.0000 0.0000 0.0000 1.3003 0.0000 13.2550 11.0080 0.0349 13.4380 0.0000 1.0000 1.0000 !OThr 1/8 epsilon i 0.0354 4.5650 0.4668 0.4668 0.8709 0.8709 0.0000 0.0000 0.0000 0.0000 6.2755 6.2755 0.0232 0.0232 1 0.0000 0.0000 0.0000 0.0000 2.0870 0.0000 0.0000 0.0000 0.0538 0.0538 0.0000 5.0000 1.7567 1.7567 0.0000 0.0000 0.0000 0.0000 0.0000 7.6397 0.1745 1.1729 13.6680 0.0000 1.0000 1.0000 0.1567 4.6094 0.4496 0.1023 -0.9998 -0.1723 0.0000 0.0000 0.0000 0.0000 2.4369 3.9490 -0.8774 0.3148 1 0.0000 0.0000 0.0000 0.0000 2.8745 0.0000 0.0000 0.0000 3.0813 0.2333 0.0000 5.0000 1.8847 2.0403 0.0000 0.0000 0.0000 1.7847 0.0000 12.0140 7.4670 6.1026 13.6260 0.0000 1.0000 1.0000 0.6738 4.5858 -0.0217 0.3581 -0.9436 0.8502 0.0000 0.0000 0.0000 0.0000 3.1814 5.1114 -0.9830 -0.9396 1 0.0000 0.0000 0.0000 0.0000 3.0907 0.0000 0.0000 0.0000 1.2890 0.1120 0.0000 5.0000 0.4529 1.6272 0.0000 0.0000 0.0000 1.6273 0.0000 7.4321 0.0349 1.5069 13.7270 0.0000 1.0000 1.0000 diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index fa701a0..e4d0f76 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -3159,6 +3159,7 @@ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 & .or. itype(i+3,1).eq.ntyp1 & .or. itype(i+4,1).eq.ntyp1) cycle +! print *,"before2",i,i+3, gshieldc_t4(2,i+3),gshieldc_t4(2,i) dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3201,6 +3202,7 @@ call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) & call eturn4(i,eello_turn4) +! print *,"before",i,i+3, gshieldc_t4(2,i+3),gshieldc_t4(2,i) num_cont_hb(i)=num_conti enddo ! i ! @@ -4545,8 +4547,10 @@ integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,& rlocshield - + j=i+3 +! if (j.ne.20) return +! print *,i,j,gshieldc_t4(2,j),gshieldc_t4(2,j+1) !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! ! Fourth-order contributions @@ -4625,6 +4629,7 @@ iresshield=shield_list(ilist,i) do k=1,3 rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i) +! print *,"rlocshield",rlocshield,grad_shield_side(k,ilist,i),iresshield gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ & rlocshield & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i) @@ -4635,16 +4640,17 @@ do ilist=1,ishield_list(j) iresshield=shield_list(ilist,j) do k=1,3 +! print *,"rlocshieldj",j,rlocshield,grad_shield_side(k,ilist,j),iresshield rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j) 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 +! print *,"after", gshieldc_t4(k,iresshield-1),iresshield-1,gshieldc_t4(k,iresshield) enddo enddo - do k=1,3 gshieldc_t4(k,i)=gshieldc_t4(k,i)+ & grad_shield(k,i)*eello_t4/fac_shield(i) @@ -4654,6 +4660,7 @@ 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) +! print *,"gshieldc_t4(k,j+1)",j,gshieldc_t4(k,j+1) enddo endif @@ -4772,10 +4779,11 @@ 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 (j.lt.nres-1) then gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) - +! endif a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) @@ -4791,10 +4799,17 @@ call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ! write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 +! if (j.lt.nres-1) then +! print *,"juest before",j1, gcorr4_turn(l,j1) gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) - +! if (shield_mode.gt.0) then +! print *,"juest after",j1, gcorr4_turn(l,j1),gshieldc_t4(k,j1),gshieldc_loc_t4(k,j1),gel_loc_turn4(i+2) +! else +! print *,"juest after",j1, gcorr4_turn(l,j1),gel_loc_turn4(i+2) +! endif +! endif enddo gshieldc_t4(3,i)=gshieldc_t4(3,i)+ & ssgradlipi*eello_t4/4.0d0*lipscale @@ -10567,7 +10582,7 @@ wscpho*gvdwc_scpho(j,i)+ & wpeppho*gvdwc_peppho(j,i) - + @@ -10765,7 +10780,11 @@ +0.5d0*(wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)& +wvdwpsb*gvdwpsb1(j,i))& +wbond_nucl*gradb_nucl(j,i)+wsbloc*gsbloc(j,i) - +! if (i.eq.21) then +! print *,"in sum",gradc(j,i,icg),wturn4*gcorr4_turn(j,i),& +! wturn4*gshieldc_t4(j,i), & +! wturn4*gshieldc_loc_t4(j,i) +! endif ! if ((i.le.2).and.(i.ge.1)) ! print *,gradc(j,i,icg),& ! gradbufc(j,i),welec*gelc(j,i), & @@ -10871,7 +10890,8 @@ ! if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i) enddo - enddo + enddo +!#define DEBUG #ifdef DEBUG write (iout,*) "gloc before adding corr" do i=1,4*nres @@ -10893,6 +10913,7 @@ write (iout,*) i,gloc(i,icg) enddo #endif +!#undef DEBUG #ifdef MPI if (nfgtasks.gt.1) then do j=1,3 @@ -11054,7 +11075,7 @@ endif endif endif -!el#define DEBUG +!#define DEBUG #ifdef DEBUG write (iout,*) "gradc gradx gloc" do i=1,nres @@ -11062,7 +11083,7 @@ i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg) enddo #endif -!el#undef DEBUG +!#undef DEBUG #ifdef TIMING time_sumgradient=time_sumgradient+MPI_Wtime()-time01 #endif @@ -11744,6 +11765,7 @@ c(j,k)=c(j,k)+aincr c(j,k+nres)=c(j,k+nres)+aincr enddo + call zerograd call etotal(energia1) etot1=energia1(0) ggg(j)=(etot1-etot)/aincr @@ -11756,6 +11778,7 @@ do j=1,3 c(j,i+nres)=c(j,i+nres)+aincr dc(j,i+nres)=dc(j,i+nres)+aincr + call zerograd call etotal(energia1) etot1=energia1(0) ggg(j+3)=(etot1-etot)/aincr @@ -11816,6 +11839,7 @@ write (iout,*) "split_ene ",split_ene call flush(iout) if (.not.split_ene) then + call zerograd call etotal(energia) etot=energia(0) call cartgrad @@ -11892,6 +11916,7 @@ dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot1=energia1(0) write (iout,*) "ij",i,j," etot1",etot1 @@ -11912,6 +11937,7 @@ dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot2=energia1(0) write (iout,*) "ij",i,j," etot2",etot2 @@ -11943,6 +11969,7 @@ dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot1=energia1(0) else @@ -11957,7 +11984,8 @@ dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then - call etotal(energia1) + call zerograd + call etotal(energia1) etot2=energia1(0) ggg(j+3)=(etot1-etot2)/(2*aincr) else @@ -12050,6 +12078,8 @@ do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) +! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i) + ! if (i.le.2) print *,"tu?!",gcart(j,i),grad_s(j,i),gxcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo @@ -12071,6 +12101,7 @@ do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) +! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo enddo @@ -12114,6 +12145,7 @@ #endif ! call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot1=energia1(0) ! call enerprint(energia1) @@ -12131,6 +12163,7 @@ call chainbuild_cart ! call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot2=energia1(0) ggg(j)=(etot1-etot2)/(2*aincr) @@ -12161,6 +12194,7 @@ ! & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2) ! write (iout,*) if (.not.split_ene) then + call zerograd call etotal(energia1) etot1=energia1(0) else @@ -12183,6 +12217,7 @@ ! write (iout,*) "dxnormnormsafe",dsqrt( ! & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2) if (.not.split_ene) then + call zerograd call etotal(energia1) etot2=energia1(0) ggg(j+3)=(etot1-etot2)/(2*aincr) @@ -16212,6 +16247,7 @@ call sum_gradient #ifdef TIMING #endif +!#define DEBUG !el write (iout,*) "After sum_gradient" #ifdef DEBUG !el write (iout,*) "After sum_gradient" @@ -16220,6 +16256,7 @@ write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) enddo #endif +!#undef DEBUG ! If performing constraint dynamics, add the gradients of the constraint energy if(usampl.and.totT.gt.eq_time) then do i=1,nct @@ -16246,6 +16283,7 @@ #endif ! call checkintcartgrad ! write(iout,*) 'calling int_to_cart' +!#define DEBUG #ifdef DEBUG write (iout,*) "gcart, gxcart, gloc before int_to_cart" #endif @@ -16277,6 +16315,7 @@ (gxcart(j,i),j=1,3) enddo #endif +!#undef DEBUG #ifdef CARGRAD #ifdef DEBUG write (iout,*) "CARGRAD" @@ -16601,22 +16640,22 @@ dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+ & cost1*dc_norm(j,i-2))/ & vbld(i-1) - domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i) + domicron(j,1,1,i)=-1.0/sint1*dcosomicron(j,1,1,i) dcosomicron(j,1,2,i)=-(dc_norm(j,i-2) & +cost1*(dc_norm(j,i-1+nres)))/ & vbld(i-1+nres) - domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i) + domicron(j,1,2,i)=-1.0/sint1*dcosomicron(j,1,2,i) !C Calculate derivative over second omicron Sci-1,Cai-1 Cai !C Looks messy but better than if in loop dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres) & +cost2*dc_norm(j,i-1))/ & vbld(i) - domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i) + domicron(j,2,1,i)=-1.0/sint2*dcosomicron(j,2,1,i) dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) & +cost2*(-dc_norm(j,i-1+nres)))/ & vbld(i-1+nres) ! write(iout,*) "vbld", i,itype(i,1),vbld(i-1+nres) - domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i) + domicron(j,2,2,i)=-1.0/sint2*dcosomicron(j,2,2,i) enddo endif enddo @@ -16677,15 +16716,20 @@ dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* & dc_norm(j,i-3))/vbld(i-2) - dphi(j,1,i)=-1/sing*dcosphi(j,1,i) + dphi(j,1,i)=-1.0/sing*dcosphi(j,1,i) dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2* & dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4* & dcostheta(j,1,i) - dphi(j,2,i)=-1/sing*dcosphi(j,2,i) + dphi(j,2,i)=-1.0/sing*dcosphi(j,2,i) dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4* & dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp* & dc_norm(j,i-1))/vbld(i) - dphi(j,3,i)=-1/sing*dcosphi(j,3,i) + dphi(j,3,i)=-1.0/sing*dcosphi(j,3,i) +!#define DEBUG +#ifdef DEBUG + write(iout,*) "just after",dphi(j,3,i),sing,dcosphi(j,3,i) +#endif +!#undef DEBUG endif enddo endif @@ -17023,6 +17067,7 @@ call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),& MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,& king,FG_COMM,IERROR) +!#define DEBUG #ifdef DEBUG !d write (iout,*) "Gather dphi" !d call flush(iout) @@ -17031,6 +17076,7 @@ write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3) enddo #endif +!#undef DEBUG call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),& MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,& king,FG_COMM,IERROR) @@ -17048,6 +17094,7 @@ #endif endif #endif +!#define DEBUG #ifdef DEBUG write (iout,*) "dtheta after gather" do i=1,nres @@ -17066,6 +17113,7 @@ write (iout,'(i3,3(3f8.5,3x))') i,((domega(j,k,i),j=1,3),k=1,3) enddo #endif +!#undef DEBUG return end subroutine intcartderiv !----------------------------------------------------------------------------- @@ -19367,7 +19415,7 @@ enddo !C now sscale fraction sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield -!C print *,buff_shield,"buff" +! print *,buff_shield,"buff",sh_frac_dist !C now sscale if (sh_frac_dist.le.0.0) cycle !C print *,ishield_list(i),i @@ -19401,6 +19449,7 @@ long=long_r_sidechain(itype(k,1)) costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2) sinthet=short/dist_pep_side*costhet +! print *,"SORT",short,long,sinthet,costhet !C now costhet_grad !C costhet=0.6d0 !C sinthet=0.8 @@ -19450,7 +19499,7 @@ enddo !C print *,sinphi,sinthet VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet)) & - & /VSolvSphere_div + /VSolvSphere_div !C & *wshield !C now the gradient... do j=1,3 @@ -19472,19 +19521,22 @@ sinphi/sinthet*costhet*costhet_grad(j)& +sinthet/sinphi*cosphi*cosphi_grad_long(j))) & )*wshield - +! print *, 1.0d0/(-dsqrt(1.0d0-sinphi*sinthet)),& +! sinphi/sinthet,& +! +sinthet/sinphi,"HERE" 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 +! print *,grad_shield_loc(j,ishield_list(i),i) enddo VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist enddo fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield) -!C write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i) +! write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i) enddo return end subroutine set_shield_fac2 diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index 2287f0a..f41bd30 100644 --- a/source/unres/geometry.f90 +++ b/source/unres/geometry.f90 @@ -3331,13 +3331,38 @@ enddo endif ! Settind dE/ddnres-1 +!#define DEBUG +#ifdef DEBUG + j=1 + write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), & + gloc(2*nres-5,icg),dtheta(j,2,nres) + +#endif +!#undef DEBUG + do j=1,3 gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ & gloc(2*nres-5,icg)*dtheta(j,2,nres) +!#define DEBUG +#ifdef DEBUG + write(iout,*)"in int to cartb",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), & + gloc(2*nres-5,icg),dtheta(j,2,nres) + +#endif +!#undef DEBUG if(itype(nres-1,1).ne.10) then gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* & dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* & domega(j,2,nres-1) +!#define DEBUG +#ifdef DEBUG + write(iout,*)"in int to cart2",i,gcart(j,nres-1),gloc(ialph(nres-1,1),icg)* & + dalpha(j,2,nres-1),gloc(ialph(nres-1,1)+nside,icg), & + domega(j,2,nres-1) + +#endif +!#undef DEBUG + endif enddo ! The side-chain vector derivatives @@ -3347,6 +3372,12 @@ do j=1,3 gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) & +gloc(ialph(i,1)+nside,icg)*domega(j,3,i) +!#define DEBUG +#ifdef DEBUG + write(iout,*)"in int to cart",i, gxcart(j,i),gloc(ialph(i,1),icg),dalpha(j,3,i), & + gloc(ialph(i,1)+nside,icg),domega(j,3,i) +#endif +!#undef DEBUG enddo endif enddo diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 7534c41..b3492fe 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -2771,6 +2771,31 @@ write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,& ' v3ss:',v3ss endif + if (shield_mode.gt.0) then + pi=4.0D0*datan(1.0D0) +!C VSolvSphere the volume of solving sphere + print *,pi,"pi" +!C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +!C there will be no distinction between proline peptide group and normal peptide +!C group in case of shielding parameters + VSolvSphere=4.0/3.0*pi*(4.50d0)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3 + write (iout,*) VSolvSphere,VSolvSphere_div +!C long axis of side chain + do i=1,ntyp + long_r_sidechain(i)=vbldsc0(1,i) +! if (scelemode.eq.0) then + short_r_sidechain(i)=sigma(i,i)/sqrt(2.0) + if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2 +! else +! short_r_sidechain(i)=sigma(i,i) +! endif + write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& + sigma0(i) + enddo + buff_shield=1.0d0 + endif + return 111 write (iout,*) "Error reading bending energy parameters." goto 999 @@ -3891,22 +3916,22 @@ write (iout,'(2a)') diagmeth(kdiag),& ' routine used to diagonalize matrices.' if (shield_mode.gt.0) then - pi=3.141592d0 + pi=4.0D0*datan(1.0D0) !C VSolvSphere the volume of solving sphere -!C print *,pi,"pi" + print *,pi,"pi" !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it !C there will be no distinction between proline peptide group and normal peptide !C group in case of shielding parameters - VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 - VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 + VSolvSphere=4.0/3.0*pi*(4.50d0)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3 write (iout,*) VSolvSphere,VSolvSphere_div !C long axis of side chain - do i=1,ntyp - long_r_sidechain(i)=vbldsc0(1,i) - short_r_sidechain(i)=sigma0(i) - write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& - sigma0(i) - enddo +! do i=1,ntyp +! long_r_sidechain(i)=vbldsc0(1,i) +! short_r_sidechain(i)=sigma0(i) +! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& +! sigma0(i) +! enddo buff_shield=1.0d0 endif return