X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fparmread.F;h=750f03e1fd0e37545817e82f91da0fbaf3c09120;hb=a1081b395546875bbbdb59cdd25f52d974adf2a6;hp=bae0f0920b9d46022e988863d0317fe48829b193;hpb=bec8e01821360e63bb9482c6c97c18fee3e97dc4;p=unres.git diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F index bae0f09..750f03e 100644 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@ -35,6 +35,7 @@ C character*16 key integer iparm double precision ip,mp +C write (iout,*) "KURWA" C C Body C @@ -55,6 +56,61 @@ C Assign virtual-bond length write (iout,*) "iparm",iparm," myparm",myparm c If reading not own parameters, skip assignment + call reada(weightcard,"D0CM",d0cm,3.78d0) + call reada(weightcard,"AKCM",akcm,15.1d0) + call reada(weightcard,"AKTH",akth,11.0d0) + call reada(weightcard,"AKCT",akct,12.0d0) + call reada(weightcard,"V1SS",v1ss,-1.08d0) + call reada(weightcard,"V2SS",v2ss,7.61d0) + call reada(weightcard,"V3SS",v3ss,13.7d0) + call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(controlcard,"DTRISS",dtriss,1.0D0) + call reada(controlcard,"ATRISS",atriss,0.3D0) + call reada(controlcard,"BTRISS",btriss,0.02D0) + call reada(controlcard,"CTRISS",ctriss,1.0D0) + dyn_ss=(index(controlcard,'DYN_SS').gt.0) +C do i=1,maxres +C dyn_ss_mask(i)=.false. +C enddo +C ebr=-12.0D0 +c +c Old arbitrary potential - commented out. +c +c dbr= 4.20D0 +c fbr= 3.30D0 +c +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 + + do i=1,maxres-1 + do j=i+1,maxres + dyn_ssbond_ij(i,j)=1.0d300 + enddo + enddo + call reada(controlcard,"HT",Ht,0.0D0) +C if (dyn_ss) then +C ss_depth=ebr/wsc-0.25*eps(1,1) +C write(iout,*) HT,wsc,eps(1,1),'KURWA' +C Ht=Ht/wsc-0.25*eps(1,1) + +C akcm=akcm*whpb/wsc +C akth=akth*whpb/wsc +C akct=akct*whpb/wsc +C v1ss=v1ss*whpb/wsc +C v2ss=v2ss*whpb/wsc +C v3ss=v3ss*whpb/wsc +C else +C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb +C endif if (iparm.eq.myparm .or. .not.separate_parset) then @@ -78,7 +134,7 @@ c wvdwpp=ww(16) wbond=ww(18) wsccor=ww(19) - + whpb=ww(15) endif call card_concat(controlcard,.false.) @@ -102,7 +158,7 @@ c Return if not own parameters call reads(controlcard,"TORDPAR",tordname_t,tordname) open (itordp,file=tordname_t,status='old') rewind(itordp) - call reads(controlcard,"SCCORAR",sccorname_t,sccorname) + call reads(controlcard,"SCCORPAR",sccorname_t,sccorname) open (isccor,file=sccorname_t,status='old') rewind(isccor) call reads(controlcard,"FOURIER",fouriername_t,fouriername) @@ -295,50 +351,61 @@ C C Read the parameters of Utheta determined from ab initio surfaces C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 C +c write (iout,*) "tu dochodze" read (ithep,*) nthetyp,ntheterm,ntheterm2, & ntheterm3,nsingle,ndouble nntheterm=max0(ntheterm,ntheterm2,ntheterm3) read (ithep,*) (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 +c write (iout,*) "tu dochodze" + 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)') res1,res2,res3 - read (ithep,*) aa0thet(i,j,k) - read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm) + enddo + do iblock=1,2 + do i=0,nthetyp + do j=-nthetyp,nthetyp + do k=-nthetyp,nthetyp + read (ithep,'(6a)') res1 + read (ithep,*) aa0thet(i,j,k,iblock) + read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm) read (ithep,*) - & ((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,*) - & (((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 @@ -350,17 +417,54 @@ C 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 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 @@ -372,16 +476,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 @@ -390,8 +494,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 @@ -690,28 +796,104 @@ C Martix of D parameters for two dimesional fourier series enddo endif #endif +C Read of Side-chain backbone correlation parameters +C Modified 11 May 2012 by Adasko +CCC C -C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local -C correlation energies. -C - read (isccor,*) nterm_sccor - do i=1,20 - do j=1,20 - read (isccor,'(a)') - do k=1,nterm_sccor - read (isccor,*) - & kk,v1sccor(k,i,j),v2sccor(k,i,j) + read (isccor,*) nsccortyp + read (isccor,*) (isccortyp(i),i=1,ntyp) + do i=-ntyp,-1 + isccortyp(i)=-isccortyp(-i) + enddo + iscprol=isccortyp(20) +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,*) + &nterm_sccor(i,j),nlor_sccor(i,j) + write (iout,*) nterm_sccor(i,j) + v0ijsccor=0.0d0 + v0ijsccor1=0.0d0 + v0ijsccor2=0.0d0 + v0ijsccor3=0.0d0 + si=-1.0d0 + nterm_sccor(-i,j)=nterm_sccor(i,j) + nterm_sccor(-i,-j)=nterm_sccor(i,j) + nterm_sccor(i,-j)=nterm_sccor(i,j) + write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j), + & nterm_sccor(-i,-j),nterm_sccor(i,-j) + do k=1,nterm_sccor(i,j) + read (isccor,*) kk,v1sccor(k,l,i,j) + & ,v2sccor(k,l,i,j) + if (j.eq.iscprol) then + if (i.eq.isccortyp(10)) then + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + else + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 + & +v2sccor(k,l,i,j)*dsqrt(0.75d0) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 + & +v1sccor(k,l,i,j)*dsqrt(0.75d0) + v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) + endif + else + if (i.eq.isccortyp(10)) then + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + else + if (j.eq.isccortyp(10)) then + v1sccor(k,l,-i,j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j) + else + v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) + endif + endif + endif + v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) + v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j) + v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j) + v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j) + si=-si + enddo + do k=1,nlor_sccor(i,j) + read (isccor,*) 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(l,i,j)=v0ijsccor + v0sccor(l,-i,j)=v0ijsccor1 + v0sccor(l,i,-j)=v0ijsccor2 + v0sccor(l,-i,-j)=v0ijsccor3 enddo enddo enddo close (isccor) if (lprint) then write (iout,'(/a/)') 'Torsional constants of SCCORR:' - do i=1,20 - do j=1,20 + do i=1,nsccortyp + do j=1,nsccortyp write (iout,*) 'ityp',i,' jtyp',j - do k=1,nterm_sccor - write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j) + write (iout,*) 'Fourier constants' + do k=1,nterm_sccor(i,j) + write (iout,'(2(1pe15.5))') + & v1sccor(k,l,i,j),v2sccor(k,l,i,j) + enddo + write (iout,*) 'Lorenz constants' + do k=1,nlor_sccor(i,j) + write (iout,'(3(1pe15.5))') + & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j) enddo enddo enddo @@ -856,8 +1038,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 + lprint=.true. if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j), & ael6(i,j),ael3(i,j) + lprint=.false. enddo enddo C @@ -898,13 +1082,17 @@ C----------------------- LJK potential -------------------------------- endif goto 50 C---------------------- GB or BP potential ----------------------------- - 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp), - & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp), - & (alp(i),i=1,ntyp) + 30 do i=1,ntyp + read (isidep,*)(eps(i,j),j=i,ntyp) + enddo + read (isidep,*)(sigma0(i),i=1,ntyp) + read (isidep,*)(sigii(i),i=1,ntyp) + read (isidep,*)(chip(i),i=1,ntyp) + read (isidep,*)(alp(i),i=1,ntyp) C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp - chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0) + chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0) enddo endif if (lprint) then @@ -1059,7 +1247,7 @@ C C C Define the constants of the disulfide bridge C - ebr=-5.50D0 +C ebr=-12.0D0 c c Old arbitrary potential - commented out. c @@ -1070,21 +1258,36 @@ 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 - if (lprint) then + if (dyn_ss) then + ss_depth=ebr/wsc-0.25*eps(1,1) +C write(iout,*) akcm,whpb,wsc,'KURWA' + Ht=Ht/wsc-0.25*eps(1,1) + + akcm=akcm*whpb/wsc + akth=akth*whpb/wsc + akct=akct*whpb/wsc + v1ss=v1ss*whpb/wsc + v2ss=v2ss*whpb/wsc + v3ss=v3ss*whpb/wsc + else + ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb + endif + +C if (lprint) 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 endif return end