X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fparmread.F;h=a38147fee5765193c998e045a5135553fdf4d7c6;hb=c1827efa66e69c93c6e2b2e7420b06b430c3550a;hp=1102c8966bec88ddde2e02559c7797043383a0f0;hpb=d101c97dea752458d76055fdbae49c26fff03c1f;p=unres.git diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 1102c89..a38147f 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"/ @@ -99,6 +101,8 @@ c enddo endif C reading lipid parameters + write (iout,*) "iliptranpar",iliptranpar + call flush(iout) read(iliptranpar,*) pepliptran do i=1,ntyp read(iliptranpar,*) liptranene(i) @@ -217,6 +221,8 @@ C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 C read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2, & ntheterm3,nsingle,ndouble + write (iout,*) "ithep",ithep + call flush(iout) nntheterm=max0(ntheterm,ntheterm2,ntheterm3) read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1) do i=-ntyp1,-1 @@ -359,11 +365,12 @@ C Control printout of the coefficients of virtual-bond-angle potentials C if (lprint) then write (iout,'(//a)') 'Parameter of virtual-bond-angle potential' - do i=1,nthetyp+1 - do j=1,nthetyp+1 - do k=1,nthetyp+1 + do iblock=1,2 + do i=0,nthetyp + do j=-nthetyp,nthetyp + do k=-nthetyp,nthetyp write (iout,'(//4a)') - & 'Type ',onelett(i),onelett(j),onelett(k) + & 'Type ',toronelet(i),toronelet(j),toronelet(k) write (iout,'(//a,10x,a)') " l","a[l]" write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock) write (iout,'(i2,1pe15.5)') @@ -392,6 +399,7 @@ C enddo enddo enddo + enddo enddo call flush(iout) endif @@ -622,6 +630,7 @@ C Read torsional parameters C read (itorp,*,end=113,err=113) ntortyp read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp) + itortyp(ntyp1)=ntortyp do iblock=1,2 do i=-ntyp,-1 itortyp(i)=-itortyp(-i) @@ -659,8 +668,9 @@ c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) close (itorp) if (lprint) then write (iout,'(/a/)') 'Torsional constants:' - do i=1,ntortyp - do j=1,ntortyp + do iblock=1,2 + do i=0,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 write (iout,*) 'ityp',i,' jtyp',j write (iout,*) 'Fourier constants' do k=1,nterm(i,j,iblock) @@ -674,6 +684,7 @@ c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) enddo enddo enddo + enddo endif C @@ -773,6 +784,80 @@ 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=0,nkcctyp + do j=0,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) v11_chyb(l,k,j,i) + enddo + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) v21_chyb(l,k,j,i) + enddo + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) v12_chyb(l,k,j,i) + enddo + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) v22_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 + if (lprint) then +c Print valence-torsional parameters + write (iout,'(a)') + & "Parameters of the valence-torsional potentials" + do i=0,nkcctyp + do j=0,nkcctyp + write (iout,'(3a)') "Type ",toronelet(i),toronelet(j) + write (iout,'(2a20,a15)') "v_kcc","v1_chyb","v2_chyb" + do k=1,nterm_kcc(j,i) + write (iout,'(i5,f15.10,i5,2f15.10)') + & k,v1_kcc(k,j,i),1,v11_chyb(1,k,j,i),v21_chyb(1,k,j,i) + do l=2,nterm_kcc_Tb(j,i) + write (iout,'(20x,i5,2f15.10)') + & l,v11_chyb(l,k,j,i),v21_chyb(l,k,j,i) + enddo + write (iout,'(i5,f15.10,i5,2f15.10)') + & k,v2_kcc(k,j,i),1,v12_chyb(1,k,j,i),v22_chyb(1,k,j,i) + do l=2,nterm_kcc_Tb(j,i) + write (iout,'(20x,i5,2f15.10)') + & l,v12_chyb(l,k,j,i),v22_chyb(l,k,j,i) + enddo + write (iout,'(a)') + enddo + enddo + enddo + endif +C here will be the apropriate recalibrating for D-aminoacid +C read (ithetkcc,*,end=121,err=121) nkcctyp + do i=0,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 + if (lprint) then + write (iout,'(a)') + & "Parameters of the valence-only potentials" + do i=0,nkcctyp + write (iout,'(2a)') "Type ",toronelet(i) + do k=1,nbend_kcc_Tb(i) + write(iout,'(i5,f15.10)') k,v1bend_chyb(k,i) + enddo + enddo + endif C Read of Side-chain backbone correlation parameters C Modified 11 May 2012 by Adasko CCC @@ -880,7 +965,7 @@ cc maxinter is maximum interaction sites v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &(1+vlor3sccor(k,i,j)**2) enddo - v0sccor(i,j,iblock)=v0ijsccor + v0sccor(l,i,j)=v0ijsccor enddo enddo enddo @@ -889,6 +974,7 @@ cc maxinter is maximum interaction sites #endif if (lprint) then write (iout,'(/a/)') 'Torsional constants:' + do l=1,maxinter do i=1,nsccortyp do j=1,nsccortyp write (iout,*) 'ityp',i,' jtyp',j @@ -903,6 +989,7 @@ cc maxinter is maximum interaction sites enddo enddo enddo + enddo endif C @@ -914,7 +1001,29 @@ C write (iout,*) "Coefficients of the cumulants" endif read (ifourier,*) nloctyp - +#ifdef NEWCORR + read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp) + read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1) + itype2loc(ntyp1)=nloctyp + iloctyp(nloctyp)=ntyp1 +#else + do i=1,ntyp1 + itype2loc(i)=itortyp(i) + enddo + iloctyp(0)=10 + iloctyp(1)=9 + iloctyp(2)=20 + iloctyp(3)=ntyp1 +#endif + do i=1,ntyp1 + itype2loc(-i)=-itype2loc(i) + enddo + do i=1,nloctyp + iloctyp(-i)=-iloctyp(i) + enddo + write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1) + write (iout,*) "nloctyp",nloctyp, + & " iloctyp",(iloctyp(i),i=0,nloctyp) do i=0,nloctyp-1 read (ifourier,*,end=115,err=115) read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13) @@ -945,16 +1054,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) - B1tilde(2,i) =-b(5) - B1tilde(1,-i) =-b(3) - B1tilde(2,-i) =b(5) - b1tilde(1,i)=0.0d0 - b1tilde(2,i)=0.0d0 - B2(1,i) = b(2) - B2(2,i) = b(4) - B2(1,-i) =b(2) - B2(2,-i) =-b(4) +cc B1tilde(1,i) = b(3,i) +cc B1tilde(2,i) =-b(5,i) +C B1tilde(1,-i) =-b(3,i) +C B1tilde(2,-i) =b(5,i) +cc b1tilde(1,i)=0.0d0 +cc b1tilde(2,i)=0.0d0 +cc B2(1,i) = b(2,i) +cc 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 @@ -1024,7 +1133,7 @@ c ee(2,1,i)=0.0d0 c ee(1,2,i)=0.0d0 c ee(2,1,i)=ee(1,2,i) enddo -c lprint=.true. + lprint=.true. if (lprint) then do i=1,nloctyp write (iout,*) 'Type',i @@ -1046,7 +1155,7 @@ c lprint=.true. enddo enddo endif -c lprint=.false. + lprint=.false. C C Read electrostatic-interaction parameters @@ -1128,11 +1237,14 @@ C---------------------- GB or BP potential ----------------------------- C now we start reading lipid do i=1,ntyp read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp) - print *,"WARNING!!" - do j=1,ntyp - epslip(i,j)=epslip(i,j)+0.05d0 - enddo + write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp) + +C print *,"WARNING!!" +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 @@ -1248,6 +1360,24 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) endif enddo enddo + write(iout,*) "tube param" + read(itube,*) epspeptube,sigmapeptube + sigmapeptube=sigmapeptube**6 + sigeps=dsign(1.0D0,epspeptube) + epspeptube=dabs(epspeptube) + pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2 + pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube + write(iout,*) pep_aa_tube,pep_bb_tube + do i=1,ntyp + read(itube,*) epssctube,sigmasctube + sigmasctube=sigmasctube**6 + sigeps=dsign(1.0D0,epssctube) + epssctube=dabs(epssctube) + sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2 + sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube + write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i) + enddo + #ifdef OLDSCP C C Define the SC-p interaction constants (hard-coded; old style) @@ -1333,6 +1463,25 @@ 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 @@ -1354,6 +1503,8 @@ C endif 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) @@ -1404,6 +1555,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