X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fparmread.F;h=ae4d7100649e77a67a91f6cbbbd1a447c41b7bbc;hb=25618f9f83673a7063414fe1e17415d138f58da8;hp=28b4385f1459ce9596f3f05039ee87ed8a5d04cd;hpb=62e5c36dfb2d05728ddba047bae10c9cf3287e38;p=unres.git diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 28b4385..ae4d710 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -26,12 +26,14 @@ 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"/ logical lprint,LaTeX dimension blower(3,3,maxlob) - dimension b(13) +C dimension b(13) character*3 lancuch,ucase C C For printing parameters after they are read set the following in the UNRES @@ -59,7 +61,7 @@ c Read the virtual-bond parameters, masses, and moments of inertia c and Stokes' radii of the peptide group and side chains c #ifdef CRYST_BOND - read (ibond,*) vbldp0,akp,mp,ip,pstok + read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok do i=1,ntyp nbondterm(i)=1 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i) @@ -71,8 +73,9 @@ c endif enddo #else - read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok + 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) @@ -97,6 +100,12 @@ c enddo enddo endif +C reading lipid parameters + read(iliptranpar,*) pepliptran + do i=1,ntyp + read(iliptranpar,*) liptranene(i) + enddo + close(iliptranpar) #ifdef CRYST_THETA C C Read the parameters of the probability distribution/energy expression @@ -212,46 +221,65 @@ C & ntheterm3,nsingle,ndouble nntheterm=max0(ntheterm,ntheterm2,ntheterm3) read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1) - do i=1,maxthetyp - do j=1,maxthetyp - do k=1,maxthetyp - aa0thet(i,j,k)=0.0d0 + do i=-ntyp1,-1 + ithetyp(i)=-ithetyp(-i) + enddo + do iblock=1,2 + do i=-maxthetyp,maxthetyp + do j=-maxthetyp,maxthetyp + do k=-maxthetyp,maxthetyp + aa0thet(i,j,k,iblock)=0.0d0 do l=1,ntheterm - aathet(l,i,j,k)=0.0d0 + aathet(l,i,j,k,iblock)=0.0d0 enddo do l=1,ntheterm2 do m=1,nsingle - bbthet(m,l,i,j,k)=0.0d0 - ccthet(m,l,i,j,k)=0.0d0 - ddthet(m,l,i,j,k)=0.0d0 - eethet(m,l,i,j,k)=0.0d0 + bbthet(m,l,i,j,k,iblock)=0.0d0 + ccthet(m,l,i,j,k,iblock)=0.0d0 + ddthet(m,l,i,j,k,iblock)=0.0d0 + eethet(m,l,i,j,k,iblock)=0.0d0 enddo enddo do l=1,ntheterm3 do m=1,ndouble do mm=1,ndouble - ffthet(mm,m,l,i,j,k)=0.0d0 - ggthet(mm,m,l,i,j,k)=0.0d0 + ffthet(mm,m,l,i,j,k,iblock)=0.0d0 + ggthet(mm,m,l,i,j,k,iblock)=0.0d0 enddo enddo enddo enddo enddo - enddo - do i=1,nthetyp - do j=1,nthetyp - do k=1,nthetyp - read (ithep,'(3a)',end=111,err=111) res1,res2,res3 - read (ithep,*,end=111,err=111) aa0thet(i,j,k) - read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm) + enddo + enddo +c VAR:iblock means terminally blocking group 1=non-proline 2=proline + do iblock=1,2 +c VAR:ntethtyp is type of theta potentials type currently 0=glycine +c VAR:1=non-glicyne non-proline 2=proline +c VAR:negative values for D-aminoacid + do i=0,nthetyp + do j=-nthetyp,nthetyp + do k=-nthetyp,nthetyp + read (ithep,'(6a)',end=111,err=111) res1 + read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock) +c VAR: aa0thet is variable describing the average value of Foureir +c VAR: expansion series +c VAR: aathet is foureir expansion in theta/2 angle for full formula +c VAR: look at the fitting equation in Kozlowska et al., J. Phys.: +Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished + read (ithep,*,end=111,err=111) + &(aathet(l,i,j,k,iblock),l=1,ntheterm) read (ithep,*,end=111,err=111) - & ((bbthet(lll,ll,i,j,k),lll=1,nsingle), - & (ccthet(lll,ll,i,j,k),lll=1,nsingle), - & (ddthet(lll,ll,i,j,k),lll=1,nsingle), - & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2) + & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle), + & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle), + & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle), + & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle), + & ll=1,ntheterm2) read (ithep,*,end=111,err=111) - & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k), - & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k), + & (((ffthet(llll,lll,ll,i,j,k,iblock), + & ffthet(lll,llll,ll,i,j,k,iblock), + & ggthet(llll,lll,ll,i,j,k,iblock), + & ggthet(lll,llll,ll,i,j,k,iblock), & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3) enddo enddo @@ -259,21 +287,75 @@ C C C For dummy ends assign glycine-type coefficients of theta-only terms; the C coefficients of theta-and-gamma-dependent terms are zero. -C +C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT +C RECOMENTDED AFTER VERSION 3.3) +c do i=1,nthetyp +c do j=1,nthetyp +c do l=1,ntheterm +c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock) +c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock) +c enddo +c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock) +c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock) +c enddo +c do l=1,ntheterm +c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock) +c enddo +c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock) +c enddo +c enddo +C AND COMMENT THE LOOPS BELOW do i=1,nthetyp do j=1,nthetyp do l=1,ntheterm - aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1) - aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j) + aathet(l,i,j,nthetyp+1,iblock)=0.0d0 + aathet(l,nthetyp+1,i,j,iblock)=0.0d0 enddo - aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1) - aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j) + aa0thet(i,j,nthetyp+1,iblock)=0.0d0 + aa0thet(nthetyp+1,i,j,iblock)=0.0d0 enddo do l=1,ntheterm - aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1) + aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0 enddo - aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1) + aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0 enddo + enddo +C TILL HERE +C Substitution for D aminoacids from symmetry. + do iblock=1,2 + do i=-nthetyp,0 + do j=-nthetyp,nthetyp + do k=-nthetyp,nthetyp + aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock) + do l=1,ntheterm + aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock) + enddo + do ll=1,ntheterm2 + do lll=1,nsingle + bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock) + ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock) + ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock) + eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock) + enddo + enddo + do ll=1,ntheterm3 + do lll=2,ndouble + do llll=1,lll-1 + ffthet(llll,lll,ll,i,j,k,iblock)= + & ffthet(llll,lll,ll,-i,-j,-k,iblock) + ffthet(lll,llll,ll,i,j,k,iblock)= + & ffthet(lll,llll,ll,-i,-j,-k,iblock) + ggthet(llll,lll,ll,i,j,k,iblock)= + & -ggthet(llll,lll,ll,-i,-j,-k,iblock) + ggthet(lll,llll,ll,i,j,k,iblock)= + & -ggthet(lll,llll,ll,-i,-j,-k,iblock) + enddo !ll + enddo !lll + enddo !llll + enddo !k + enddo !j + enddo !i + enddo !iblock C C Control printout of the coefficients of virtual-bond-angle potentials C @@ -285,16 +367,16 @@ C write (iout,'(//4a)') & 'Type ',onelett(i),onelett(j),onelett(k) write (iout,'(//a,10x,a)') " l","a[l]" - write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k) + write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock) write (iout,'(i2,1pe15.5)') - & (l,aathet(l,i,j,k),l=1,ntheterm) + & (l,aathet(l,i,j,k,iblock),l=1,ntheterm) do l=1,ntheterm2 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') & "b",l,"c",l,"d",l,"e",l do m=1,nsingle write (iout,'(i2,4(1pe15.5))') m, - & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k), - & ddthet(m,l,i,j,k),eethet(m,l,i,j,k) + & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock), + & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock) enddo enddo do l=1,ntheterm3 @@ -303,8 +385,10 @@ C do m=2,ndouble do n=1,m-1 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m, - & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k), - & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k) + & ffthet(n,m,l,i,j,k,iblock), + & ffthet(m,n,l,i,j,k,iblock), + & ggthet(n,m,l,i,j,k,iblock), + & ggthet(m,n,l,i,j,k,iblock) enddo enddo enddo @@ -313,6 +397,53 @@ C enddo call flush(iout) endif + write (2,*) "Start reading THETA_PDB",ithep_pdb + do i=1,ntyp +c write (2,*) 'i=',i + read (ithep_pdb,*,err=111,end=111) + & a0thet(i),(athet(j,i,1,1),j=1,2), + & (bthet(j,i,1,1),j=1,2) + read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3) + read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3) + read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i) + sigc0(i)=sigc0(i)**2 + enddo + do i=1,ntyp + athet(1,i,1,-1)=athet(1,i,1,1) + athet(2,i,1,-1)=athet(2,i,1,1) + bthet(1,i,1,-1)=-bthet(1,i,1,1) + bthet(2,i,1,-1)=-bthet(2,i,1,1) + athet(1,i,-1,1)=-athet(1,i,1,1) + athet(2,i,-1,1)=-athet(2,i,1,1) + bthet(1,i,-1,1)=bthet(1,i,1,1) + bthet(2,i,-1,1)=bthet(2,i,1,1) + enddo + do i=-ntyp,-1 + a0thet(i)=a0thet(-i) + athet(1,i,-1,-1)=athet(1,-i,1,1) + athet(2,i,-1,-1)=-athet(2,-i,1,1) + bthet(1,i,-1,-1)=bthet(1,-i,1,1) + bthet(2,i,-1,-1)=-bthet(2,-i,1,1) + athet(1,i,-1,1)=athet(1,-i,1,1) + athet(2,i,-1,1)=-athet(2,-i,1,1) + bthet(1,i,-1,1)=-bthet(1,-i,1,1) + bthet(2,i,-1,1)=bthet(2,-i,1,1) + athet(1,i,1,-1)=-athet(1,-i,1,1) + athet(2,i,1,-1)=athet(2,-i,1,1) + bthet(1,i,1,-1)=bthet(1,-i,1,1) + bthet(2,i,1,-1)=-bthet(2,-i,1,1) + theta0(i)=theta0(-i) + sig0(i)=sig0(-i) + sigc0(i)=sigc0(-i) + do j=0,3 + polthet(j,i)=polthet(j,-i) + enddo + do j=1,3 + gthet(j,i)=gthet(j,-i) + enddo + enddo + write (2,*) "End reading THETA_PDB" + close (ithep_pdb) #endif close(ithep) #ifdef CRYST_SC @@ -415,6 +546,50 @@ C enddo endif enddo +C +C Read the parameters of the probability distribution/energy expression +C of the side chains. +C + write (2,*) "Start reading ROTAM_PDB" + do i=1,ntyp + read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i) + if (i.eq.10) then + dsc_inv(i)=0.0D0 + else + dsc_inv(i)=1.0D0/dsc(i) + endif + if (i.ne.10) then + do j=1,nlob(i) + do k=1,3 + do l=1,3 + blower(l,k,j)=0.0D0 + enddo + enddo + enddo + bsc(1,i)=0.0D0 + read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3), + & ((blower(k,l,1),l=1,k),k=1,3) + do j=2,nlob(i) + read (irotam_pdb,*,end=112,err=112) bsc(j,i) + read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3), + & ((blower(k,l,j),l=1,k),k=1,3) + enddo + do j=1,nlob(i) + do k=1,3 + do l=1,k + akl=0.0D0 + do m=1,3 + akl=akl+blower(k,m,j)*blower(l,m,j) + enddo + gaussc(k,l,j,i)=akl + gaussc(l,k,j,i)=akl + enddo + enddo + enddo + endif + enddo + close (irotam_pdb) + write (2,*) "End reading ROTAM_PDB" #endif close(irotam) @@ -669,7 +844,7 @@ cc maxinter is maximum interaction sites si=-si enddo do k=1,nlor_sccor(i,j) - read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j), + read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j), & vlor2sccor(k,i,j),vlor3sccor(k,i,j) v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &(1+vlor3sccor(k,i,j)**2) @@ -683,31 +858,31 @@ cc maxinter is maximum interaction sites enddo close (isccor) #else - read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp) + read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp) c write (iout,*) 'ntortyp',ntortyp maxinter=3 cc maxinter is maximum interaction sites do l=1,maxinter do i=1,nsccortyp do j=1,nsccortyp - read (isccor,*,end=113,err=113) + read (isccor,*,end=119,err=119) & nterm_sccor(i,j),nlor_sccor(i,j) v0ijsccor=0.0d0 si=-1.0d0 do k=1,nterm_sccor(i,j) - read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j) + read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j) & ,v2sccor(k,l,i,j) v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) si=-si enddo do k=1,nlor_sccor(i,j) - read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j), + read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j), & vlor2sccor(k,i,j),vlor3sccor(k,i,j) v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &(1+vlor3sccor(k,i,j)**2) enddo - v0sccor(i,j)=v0ijsccor + v0sccor(i,j,iblock)=v0ijsccor enddo enddo enddo @@ -741,97 +916,117 @@ C write (iout,*) "Coefficients of the cumulants" endif read (ifourier,*) nloctyp + do i=0,nloctyp-1 read (ifourier,*,end=115,err=115) - read (ifourier,*,end=115,err=115) (b(ii),ii=1,13) + read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13) +#ifdef NEWCORR + read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3) + read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3) + read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1) + read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1) + read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1) +#endif if (lprint) then write (iout,*) 'Type',i - write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13) + write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13) endif - B1(1,i) = b(3) - B1(2,i) = b(5) - B1(1,-i) = b(3) - B1(2,-i) = -b(5) +c B1(1,i) = b(3) +c B1(2,i) = b(5) +c B1(1,-i) = b(3) +c B1(2,-i) = -b(5) c b1(1,i)=0.0d0 c b1(2,i)=0.0d0 - B1tilde(1,i) = b(3) - B1tilde(2,i) =-b(5) - B1tilde(1,-i) =-b(3) - B1tilde(2,-i) =b(5) +c B1tilde(1,i) = b(3) +c B1tilde(2,i) =-b(5) +c B1tilde(1,-i) =-b(3) +c B1tilde(2,-i) =b(5) c b1tilde(1,i)=0.0d0 c 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) +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 - CC(1,1,i)= b(7) - CC(2,2,i)=-b(7) - CC(2,1,i)= b(9) - CC(1,2,i)= b(9) - CC(1,1,-i)= b(7) - CC(2,2,-i)=-b(7) - CC(2,1,-i)=-b(9) - CC(1,2,-i)=-b(9) + CC(1,1,i)= b(7,i) + CC(2,2,i)=-b(7,i) + CC(2,1,i)= b(9,i) + CC(1,2,i)= b(9,i) + CC(1,1,-i)= b(7,i) + CC(2,2,-i)=-b(7,i) + CC(2,1,-i)=-b(9,i) + CC(1,2,-i)=-b(9,i) c CC(1,1,i)=0.0d0 c CC(2,2,i)=0.0d0 c CC(2,1,i)=0.0d0 c CC(1,2,i)=0.0d0 - Ctilde(1,1,i)=b(7) - Ctilde(1,2,i)=b(9) - Ctilde(2,1,i)=-b(9) - Ctilde(2,2,i)=b(7) - Ctilde(1,1,-i)=b(7) - Ctilde(1,2,-i)=-b(9) - Ctilde(2,1,-i)=b(9) - Ctilde(2,2,-i)=b(7) + Ctilde(1,1,i)=b(7,i) + Ctilde(1,2,i)=b(9,i) + Ctilde(2,1,i)=-b(9,i) + Ctilde(2,2,i)=b(7,i) + Ctilde(1,1,-i)=b(7,i) + Ctilde(1,2,-i)=-b(9,i) + Ctilde(2,1,-i)=b(9,i) + Ctilde(2,2,-i)=b(7,i) c Ctilde(1,1,i)=0.0d0 c Ctilde(1,2,i)=0.0d0 c Ctilde(2,1,i)=0.0d0 c Ctilde(2,2,i)=0.0d0 - DD(1,1,i)= b(6) - DD(2,2,i)=-b(6) - DD(2,1,i)= b(8) - DD(1,2,i)= b(8) - DD(1,1,-i)= b(6) - DD(2,2,-i)=-b(6) - DD(2,1,-i)=-b(8) - DD(1,2,-i)=-b(8) + DD(1,1,i)= b(6,i) + DD(2,2,i)=-b(6,i) + DD(2,1,i)= b(8,i) + DD(1,2,i)= b(8,i) + DD(1,1,-i)= b(6,i) + DD(2,2,-i)=-b(6,i) + DD(2,1,-i)=-b(8,i) + DD(1,2,-i)=-b(8,i) c DD(1,1,i)=0.0d0 c DD(2,2,i)=0.0d0 c DD(2,1,i)=0.0d0 c DD(1,2,i)=0.0d0 - Dtilde(1,1,i)=b(6) - Dtilde(1,2,i)=b(8) - Dtilde(2,1,i)=-b(8) - Dtilde(2,2,i)=b(6) - Dtilde(1,1,-i)=b(6) - Dtilde(1,2,-i)=-b(8) - Dtilde(2,1,-i)=b(8) - Dtilde(2,2,-i)=b(6) + Dtilde(1,1,i)=b(6,i) + Dtilde(1,2,i)=b(8,i) + Dtilde(2,1,i)=-b(8,i) + Dtilde(2,2,i)=b(6,i) + Dtilde(1,1,-i)=b(6,i) + Dtilde(1,2,-i)=-b(8,i) + Dtilde(2,1,-i)=b(8,i) + Dtilde(2,2,-i)=b(6,i) c Dtilde(1,1,i)=0.0d0 c Dtilde(1,2,i)=0.0d0 c Dtilde(2,1,i)=0.0d0 c Dtilde(2,2,i)=0.0d0 - EE(1,1,i)= b(10)+b(11) - EE(2,2,i)=-b(10)+b(11) - EE(2,1,i)= b(12)-b(13) - EE(1,2,i)= b(12)+b(13) - EE(1,1,-i)= b(10)+b(11) - EE(2,2,-i)=-b(10)+b(11) - EE(2,1,-i)=-b(12)+b(13) - EE(1,2,-i)=-b(12)-b(13) - + EEold(1,1,i)= b(10,i)+b(11,i) + EEold(2,2,i)=-b(10,i)+b(11,i) + EEold(2,1,i)= b(12,i)-b(13,i) + EEold(1,2,i)= b(12,i)+b(13,i) + EEold(1,1,-i)= b(10,i)+b(11,i) + EEold(2,2,-i)=-b(10,i)+b(11,i) + 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 c ee(1,2,i)=0.0d0 c ee(2,1,i)=ee(1,2,i) enddo +c lprint=.true. if (lprint) then do i=1,nloctyp write (iout,*) 'Type',i @@ -849,10 +1044,11 @@ c ee(2,1,i)=ee(1,2,i) enddo write(iout,*) 'EE' do j=1,2 - write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i) + write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i) enddo enddo endif +c lprint=.false. C C Read electrostatic-interaction parameters @@ -875,8 +1071,10 @@ C bpp (i,j)=-2.0D0*epp(i,j)*rri ael6(i,j)=elpp6(i,j)*4.2D0**6 ael3(i,j)=elpp3(i,j)*4.2D0**3 +c lprint=.true. if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j), & ael6(i,j),ael3(i,j) +c lprint=.false. enddo enddo C @@ -897,7 +1095,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:' @@ -909,7 +1107,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:' @@ -923,12 +1121,21 @@ 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) read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp) read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp) +C now we start reading lipid + do i=1,ntyp + read (isidep,*,end=1161,err=1161)(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 @@ -947,7 +1154,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 @@ -967,6 +1174,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 @@ -995,10 +1203,17 @@ C Calculate the "working" parameters of SC interactions. epsij=eps(i,j) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) - aa(i,j)=epsij*rrij*rrij - bb(i,j)=-sigeps*epsij*rrij - aa(j,i)=aa(i,j) - bb(j,i)=bb(i,j) + aa_aq(i,j)=epsij*rrij*rrij + bb_aq(i,j)=-sigeps*epsij*rrij + aa_aq(j,i)=aa_aq(i,j) + bb_aq(j,i)=bb_aq(i,j) + epsijlip=epslip(i,j) + sigeps=dsign(1.0D0,epsijlip) + epsijlip=dabs(epsijlip) + aa_lip(i,j)=epsijlip*rrij*rrij + bb_lip(i,j)=-sigeps*epsijlip*rrij + aa_lip(j,i)=aa_lip(i,j) + bb_lip(j,i)=bb_lip(i,j) if (ipot.gt.2) then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 @@ -1031,7 +1246,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) endif if (lprint) then write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') - & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j), + & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j), & sigma(i,j),r0(i,j),chi(i,j),chi(j,i) endif enddo @@ -1075,7 +1290,7 @@ C bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6 enddo - +c lprint=.true. if (lprint) then write (iout,*) "Parameters of SC-p interactions:" do i=1,ntyp @@ -1083,11 +1298,12 @@ C & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2) enddo endif +c lprint=.false. #endif 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 @@ -1098,13 +1314,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 @@ -1112,14 +1328,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 @@ -1134,6 +1369,8 @@ c v3ss=0.0d0 goto 999 116 write (iout,*) "Error reading electrostatic energy parameters." goto 999 + 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip" + goto 999 117 write (iout,*) "Error reading side chain interaction parameters." goto 999 118 write (iout,*) "Error reading SCp interaction parameters." @@ -1189,6 +1426,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