X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fparmread.F;h=438f544935800998d7e26ca758834ae37d5d0d79;hb=88af2feafd40050318384e0a14aac77a4f9fffa4;hp=093ca8afad8acd3e0a7caf1352fde74b2dd8ed63;hpb=7bb45761e592e3631013f7f09a37c57068707c48;p=unres.git diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 093ca8a..438f544 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -26,6 +26,8 @@ C include 'COMMON.SBRIDGE' include 'COMMON.MD' include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.SHIELD' character*1 t1,t2,t3 character*1 onelett(4) /"G","A","P","D"/ character*1 toronelet(-2:2) /"p","a","G","A","P"/ @@ -73,6 +75,7 @@ c #else read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok do i=1,ntyp + print *,i read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i), & j=1,nbondterm(i)),msc(i),isc(i),restok(i) dsc(i) = vbldsc0(1,i) @@ -772,6 +775,43 @@ C Martix of D parameters for two dimesional fourier series enddo endif #endif +C read Czybyshev torsional parameters + read (itorkcc,*,end=121,err=121) nkcctyp + read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp) + do i=-ntyp,-1 + itortyp_kcc(i)=-itortyp_kcc(-i) + enddo + do i=1,nkcctyp + do j=1,nkcctyp +C first we read the cos and sin gamma parameters + read (itorkcc,*,end=121,err=121) + & nterm_kcc(j,i),nterm_kcc_Tb(j,i) +C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i) + do k=1,nterm_kcc(j,i) + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) + & v1_chyb(l,k,j,i) + enddo + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) + & v2_chyb(l,k,j,i) + enddo + read (itorkcc,*,end=121,err=121) + & v1_kcc(k,j,i) + read (itorkcc,*,end=121,err=121) + & v2_kcc(k,j,i) + enddo + enddo + enddo +C here will be the apropriate recalibrating for D-aminoacid +C read (ithetkcc,*,end=121,err=121) nkcctyp + do i=1,nkcctyp + read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i) + do j=1,nbend_kcc_Tb(i) + read (ithetkcc,*,end=121,err=121) + & v1bend_chyb(j,i) + enddo + enddo C Read of Side-chain backbone correlation parameters C Modified 11 May 2012 by Adasko CCC @@ -944,6 +984,16 @@ c B2(1,i) = b(2) c B2(2,i) = b(4) c B2(1,-i) =b(2) c B2(2,-i) =-b(4) + B1tilde(1,i) = b(3,i) + B1tilde(2,i) =-b(5,i) +C B1tilde(1,-i) =-b(3,i) +C B1tilde(2,-i) =b(5,i) + b1tilde(1,i)=0.0d0 + b1tilde(2,i)=0.0d0 + B2(1,i) = b(2,i) + B2(2,i) = b(4,i) +C B2(1,-i) =b(2,i) +C B2(2,-i) =-b(4,i) c b2(1,i)=0.0d0 c b2(2,i)=0.0d0 @@ -1082,7 +1132,7 @@ C & ', exponents are ',expon,2*expon goto (10,20,30,30,40) ipot C----------------------- LJ potential --------------------------------- - 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), + 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp), & (sigma0(i),i=1,ntyp) if (lprint) then write (iout,'(/a/)') 'Parameters of the LJ potential:' @@ -1094,7 +1144,7 @@ C----------------------- LJ potential --------------------------------- endif goto 50 C----------------------- LJK potential -------------------------------- - 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), + 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp), & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp) if (lprint) then write (iout,'(/a/)') 'Parameters of the LJK potential:' @@ -1108,7 +1158,7 @@ C----------------------- LJK potential -------------------------------- goto 50 C---------------------- GB or BP potential ----------------------------- 30 do i=1,ntyp - read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp) + read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp) enddo read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp) read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp) @@ -1122,6 +1172,7 @@ C do j=1,ntyp C epslip(i,j)=epslip(i,j)+0.05d0 C enddo enddo + write(iout,*) epslip(1,1),"OK?" C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@ -1140,7 +1191,7 @@ C For the GB potential convert sigma'**2 into chi' endif goto 50 C--------------------- GBV potential ----------------------------------- - 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), + 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp), & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp), & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp) if (lprint) then @@ -1289,7 +1340,7 @@ c lprint=.false. C C Define the constants of the disulfide bridge C - ebr=-5.50D0 +C ebr=-12.00D0 c c Old arbitrary potential - commented out. c @@ -1300,13 +1351,13 @@ c Constants of the disulfide-bond potential determined based on the RHF/6-31G** c energy surface of diethyl disulfide. c A. Liwo and U. Kozlowska, 11/24/03 c - D0CM = 3.78d0 - AKCM = 15.1d0 - AKTH = 11.0d0 - AKCT = 12.0d0 - V1SS =-1.08d0 - V2SS = 7.61d0 - V3SS = 13.7d0 +C D0CM = 3.78d0 +C AKCM = 15.1d0 +C AKTH = 11.0d0 +C AKCT = 12.0d0 +C V1SS =-1.08d0 +C V2SS = 7.61d0 +C V3SS = 13.7d0 c akcm=0.0d0 c akth=0.0d0 c akct=0.0d0 @@ -1314,14 +1365,33 @@ c v1ss=0.0d0 c v2ss=0.0d0 c v3ss=0.0d0 - if(me.eq.king) then - write (iout,'(/a)') "Disulfide bridge parameters:" - write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr - write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm - write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct - write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, - & ' v3ss:',v3ss - endif +C if(me.eq.king) then +C write (iout,'(/a)') "Disulfide bridge parameters:" +C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr +C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm +C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct +C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, +C & ' v3ss:',v3ss +C endif +C set the variables used for shielding effect +C write (iout,*) "SHIELD MODE",shield_mode +C if (shield_mode.gt.0) then +C VSolvSphere the volume of solving sphere +C 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 +C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 +C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 +C write (iout,*) VSolvSphere,VSolvSphere_div +C long axis of side chain +C do i=1,ntyp +C long_r_sidechain(i)=vbldsc0(1,i) +C short_r_sidechain(i)=sigma0(i) +C enddo +C lets set the buffor value +C buff_shield=1.0d0 +C endif return 111 write (iout,*) "Error reading bending energy parameters." goto 999 @@ -1343,6 +1413,8 @@ c v3ss=0.0d0 118 write (iout,*) "Error reading SCp interaction parameters." goto 999 119 write (iout,*) "Error reading SCCOR parameters" + go to 999 + 121 write (iout,*) "Error in Czybyshev parameters" 999 continue #ifdef MPI call MPI_Finalize(Ierror) @@ -1393,6 +1465,22 @@ c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--' #else call getenv(var,val) #endif - +C set the variables used for shielding effect +C if (shield_mode.gt.0) then +C VSolvSphere the volume of solving sphere +C 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 +C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 +C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 +C long axis of side chain +C do i=1,ntyp +C long_r_sidechain(i)=vbldsc0(1,i) +C short_r_sidechain(i)=sigma0(i) +C enddo +C lets set the buffor value +C buff_shield=1.0d0 +C endif return end