X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fparmread.F;h=2d2c23ba8d5a4bb9bbce513926b610c3b5c96569;hb=d305127baa0f7e25ecb5422e6e0e34aa15db4776;hp=9ba20446e5528ce2aba212f4f6e63f06fadde794;hpb=4ca1e88a5f4f23619f7eaa7edcc35d5c4aca9b8c;p=unres.git diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 9ba2044..2d2c23b 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) @@ -98,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) @@ -216,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 @@ -358,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)') @@ -391,6 +399,7 @@ C enddo enddo enddo + enddo enddo call flush(iout) endif @@ -658,8 +667,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) @@ -673,6 +683,7 @@ c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) enddo enddo enddo + enddo endif C @@ -772,6 +783,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 @@ -879,7 +964,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 @@ -888,6 +973,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 @@ -902,6 +988,7 @@ cc maxinter is maximum interaction sites enddo enddo enddo + enddo endif C @@ -913,7 +1000,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) @@ -944,6 +1053,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) +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 @@ -1006,6 +1125,7 @@ c Dtilde(2,2,i)=0.0d0 EEold(2,1,-i)=-b(12,i)+b(13,i) EEold(1,2,-i)=-b(12,i)-b(13,i) write(iout,*) "TU DOCHODZE" + print *,"JESTEM" c ee(1,1,i)=1.0d0 c ee(2,2,i)=1.0d0 c ee(2,1,i)=0.0d0 @@ -1081,7 +1201,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:' @@ -1093,7 +1213,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:' @@ -1107,7 +1227,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) @@ -1116,11 +1236,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.5d0 - 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 @@ -1139,7 +1262,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 @@ -1159,6 +1282,7 @@ C Calculate the "working" parameters of SC interactions. do i=2,ntyp do j=1,i-1 eps(i,j)=eps(j,i) + epslip(i,j)=epslip(j,i) enddo enddo do i=1,ntyp @@ -1287,7 +1411,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 @@ -1298,13 +1422,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 @@ -1312,14 +1436,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 @@ -1341,6 +1484,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) @@ -1391,6 +1536,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