X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fparmread.F;h=f6ea06a52205cfe6c69d0abc8fe91d5f8c570fa4;hb=18010d3bb85a86b499dc4569225e21faffea9634;hp=654c30518eb9f6067150240f4c3858f80a786fac;hpb=88ed3524aab7542a31b89909e78f9cd79b8fe99b;p=unres.git diff --git a/source/unres/src_MD/parmread.F b/source/unres/src_MD/parmread.F index 654c305..f6ea06a 100644 --- a/source/unres/src_MD/parmread.F +++ b/source/unres/src_MD/parmread.F @@ -28,6 +28,7 @@ C include 'COMMON.SETUP' 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) @@ -102,13 +103,48 @@ C Read the parameters of the probability distribution/energy expression C of the virtual-bond valence angles theta C do i=1,ntyp - read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2), - & (bthet(j,i),j=1,2) + read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2), + & (bthet(j,i,1,1),j=1,2) read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3) - read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3) - read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i) - sigc0(i)=sigc0(i)**2 + read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3) + read (ithep,*,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 + close (ithep) if (lprint) then if (.not.LaTeX) then @@ -119,7 +155,7 @@ C & ' B1 ',' B2 ' do i=1,ntyp write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i, - & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2) + & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2) enddo write (iout,'(/a/9x,5a/79(1h-))') & 'Parameters of the expression for sigma(theta_c):', @@ -146,7 +182,8 @@ C & ' b1*10^1 ',' b2*10^1 ' do i=1,ntyp write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i), - & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2) + & a0thet(i),(100*athet(j,i,1,1),j=1,2), + & (10*bthet(j,i,1,1),j=1,2) enddo write (iout,'(/a/9x,5a/79(1h-))') & 'Parameters of the expression for sigma(theta_c):', @@ -174,69 +211,146 @@ C read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2, & ntheterm3,nsingle,ndouble nntheterm=max0(ntheterm,ntheterm2,ntheterm3) + write(iout,*) "I am here",ntyp1 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 +C write(iout,*) "I am herew" + 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 +c VAR: aa0thet is variable describing the average value of Foureir +c VAR: expansion series + read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock) +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 enddo + + 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 @@ -245,19 +359,19 @@ C do i=1,nthetyp+1 do j=1,nthetyp+1 do k=1,nthetyp+1 - write (iout,'(//4a)') - & 'Type ',onelett(i),onelett(j),onelett(k) + 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]))') + 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 @@ -266,8 +380,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 @@ -278,13 +394,48 @@ C endif write (2,*) "Start reading THETA_PDB" do i=1,ntyp - read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2), - & (bthet(j,i),j=1,2) + 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 @@ -312,18 +463,10 @@ C bsc(1,i)=0.0D0 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3), & ((blower(k,l,1),l=1,k),k=1,3) - censc(1,1,-i)=censc(1,1,i) - censc(2,1,-i)=censc(2,1,i) - censc(3,1,-i)=-censc(3,1,i) - do j=2,nlob(i) read (irotam,*,end=112,err=112) bsc(j,i) read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3), & ((blower(k,l,j),l=1,k),k=1,3) - censc(1,j,-i)=censc(1,j,i) - censc(2,j,-i)=censc(2,j,i) - censc(3,j,-i)=-censc(3,j,i) -C BSC is amplitude of Gaussian enddo do j=1,nlob(i) do k=1,3 @@ -334,14 +477,6 @@ C BSC is amplitude of Gaussian enddo gaussc(k,l,j,i)=akl gaussc(l,k,j,i)=akl - if (((k.eq.3).and.(l.ne.3)) - & .or.((l.eq.3).and.(k.ne.3))) then - gaussc(k,l,j,-i)=-akl - gaussc(l,k,j,-i)=-akl - else - gaussc(k,l,j,-i)=akl - gaussc(l,k,j,-i)=akl - endif enddo enddo enddo @@ -394,6 +529,8 @@ 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 @@ -432,6 +569,7 @@ C endif enddo close (irotam_pdb) + write (2,*) "Ending reading ROTAM_PDB" #endif close(irotam) @@ -470,26 +608,29 @@ C do i=-ntyp,-1 itortyp(i)=-itortyp(-i) enddo -c write (iout,*) 'ntortyp',ntortyp + write (iout,*) 'ntortyp',ntortyp do i=0,ntortyp-1 - do j=-ntortyp,ntortyp - read (itorp,*,end=113,err=113) nterm(i,j,iblock), + do j=-ntortyp+1,ntortyp-1 + read (itorp,*,end=113,err=113) nterm(i,j,iblock), & nlor(i,j,iblock) nterm(-i,-j,iblock)=nterm(i,j,iblock) nlor(-i,-j,iblock)=nlor(i,j,iblock) v0ij=0.0d0 si=-1.0d0 - do k=1,nterm(i,j,iblock) - read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock), - & v2(k,i,j,iblock) + do k=1,nterm(i,j,iblock) + read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock), + & v2(k,i,j,iblock) v1(k,-i,-j,iblock)=v1(k,i,j,iblock) v2(k,-i,-j,iblock)=-v2(k,i,j,iblock) - v0ij=v0ij+si*v1(k,i,j) + v0ij=v0ij+si*v1(k,i,j,iblock) si=-si +c write(iout,*) i,j,k,iblock,nterm(i,j,iblock) +c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock), +c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) enddo - do k=1,nlor(i,j,iblock) + do k=1,nlor(i,j,iblock) read (itorp,*,end=113,err=113) kk,vlor1(k,i,j), - & vlor2(k,i,j),vlor3(k,i,j) + & vlor2(k,i,j),vlor3(k,i,j) v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2) enddo v0(i,j,iblock)=v0ij @@ -499,17 +640,18 @@ c write (iout,*) 'ntortyp',ntortyp enddo close (itorp) if (lprint) then - write (iout,'(/a/)') 'Torsional constants:' - do i=1,ntortyp - do j=1,ntortyp + write (iout,'(/a/)') 'Torsional constants:' + do i=1,ntortyp + do j=1,ntortyp write (iout,*) 'ityp',i,' jtyp',j write (iout,*) 'Fourier constants' - do k=1,nterm(i,j) - write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j) + do k=1,nterm(i,j,iblock) + write (iout,'(2(1pe15.5))') v1(k,i,j,iblock), + & v2(k,i,j,iblock) enddo write (iout,*) 'Lorenz constants' - do k=1,nlor(i,j) - write (iout,'(3(1pe15.5))') + do k=1,nlor(i,j,iblock) + write (iout,'(3(1pe15.5))') & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) enddo enddo @@ -518,12 +660,19 @@ c write (iout,*) 'ntortyp',ntortyp C C 6/23/01 Read parameters for double torsionals C - do i=1,ntortyp - do j=1,ntortyp - do k=1,ntortyp +C do i=1,ntortyp +C do j=1,ntortyp +C do k=1,ntortyp + do iblock=1,2 + do i=0,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 + do k=-ntortyp+1,ntortyp-1 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3 - if (t1.ne.onelett(i) .or. t2.ne.onelett(j) - & .or. t3.ne.onelett(k)) then +c write (iout,*) "OK onelett", +c & i,j,k,t1,t2,t3 + + if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) + & .or. t3.ne.toronelet(k)) then write (iout,*) "Error in double torsional parameter file", & i,j,k,t1,t2,t3 #ifdef MPI @@ -531,94 +680,202 @@ C #endif stop "Error in double torsional parameter file" endif - read (itordp,*,end=114,err=114) ntermd_1(i,j,k), - & ntermd_2(i,j,k) - read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1, - & ntermd_1(i,j,k)) - read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1, - & ntermd_1(i,j,k)) - read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1, - & ntermd_1(i,j,k)) - read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1, - & ntermd_1(i,j,k)) - read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k), - & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k), - & m=1,l-1),l=1,ntermd_2(i,j,k)) - enddo - enddo - enddo + read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock), + & ntermd_2(i,j,k,iblock) + ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock) + ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock) + read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1, + & ntermd_1(i,j,k,iblock)) + read (itordp,*,end=114,err=114)(v1s(1,l,i,j,k,iblock),l=1, + & ntermd_1(i,j,k,iblock)) + read (itordp,*,end=114,err=114)(v1c(2,l,i,j,k,iblock),l=1, + & ntermd_1(i,j,k,iblock)) + read (itordp,*,end=114,err=114)(v1s(2,l,i,j,k,iblock),l=1, + & ntermd_1(i,j,k,iblock)) +C Martix of D parameters for one dimesional foureir series + do l=1,ntermd_1(i,j,k,iblock) + v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock) + v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock) + v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock) + v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock) +c write(iout,*) "whcodze" , +c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock) + enddo + read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock), + & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock), + & v2s(m,l,i,j,k,iblock), + & m=1,l-1),l=1,ntermd_2(i,j,k,iblock)) +C Martix of D parameters for two dimesional fourier series + do l=1,ntermd_2(i,j,k,iblock) + do m=1,l-1 + v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock) + v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock) + v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock) + v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock) + enddo!m + enddo!l + enddo!k + enddo!j + enddo!i + enddo!iblock if (lprint) then - write (iout,*) + write (iout,*) write (iout,*) 'Constants for double torsionals' - do i=1,ntortyp - do j=1,ntortyp - do k=1,ntortyp + do iblock=1,2 + do i=0,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 + do k=-ntortyp+1,ntortyp-1 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k, - & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k) + & ' nsingle',ntermd_1(i,j,k,iblock), + & ' ndouble',ntermd_2(i,j,k,iblock) write (iout,*) write (iout,*) 'Single angles:' - do l=1,ntermd_1(i,j,k) - write (iout,'(i5,2f10.5,5x,2f10.5)') l, - & v1c(1,l,i,j,k),v1s(1,l,i,j,k), - & v1c(2,l,i,j,k),v1s(2,l,i,j,k) + do l=1,ntermd_1(i,j,k,iblock) + write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l, + & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock), + & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock), + & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock) enddo write (iout,*) write (iout,*) 'Pairs of angles:' - write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k)) - do l=1,ntermd_2(i,j,k) - write (iout,'(i5,20f10.5)') - & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k)) + write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock)) + do l=1,ntermd_2(i,j,k,iblock) + write (iout,'(i5,20f10.5)') + & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)) enddo write (iout,*) - write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k)) - do l=1,ntermd_2(i,j,k) - write (iout,'(i5,20f10.5)') - & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k)) + write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock)) + do l=1,ntermd_2(i,j,k,iblock) + write (iout,'(i5,20f10.5)') + & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)), + & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock)) enddo write (iout,*) enddo enddo enddo + enddo endif #endif C Read of Side-chain backbone correlation parameters C Modified 11 May 2012 by Adasko CCC C - read (isccor,*,end=113,err=113) nsccortyp + read (isccor,*,end=1113,err=1113) nsccortyp #ifdef SCCORPDB - read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp) -#else - read (isccor,*,end=113,err=113) (isccortyp(i),i=-ntyp,ntyp) -#endif -c write (iout,*) 'ntortyp',ntortyp + write (iout,*) "Tu wchodze" + read (isccor,*,end=1113,err=1113) (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,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j) + do j=1,nsccortyp + read (isccor,*,end=1113,err=1113) nterm_sccor(i,j), + & nlor_sccor(i,j) v0ijsccor=0.0d0 + v0ijsccor1=0.0d0 + v0ijsccor2=0.0d0 + v0ijsccor3=0.0d0 si=-1.0d0 - - do k=1,nterm_sccor(i,j) - read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j) - & ,v2sccor(k,l,i,j) + nterm_sccor(-i,j)=nterm_sccor(i,j) + 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,*,end=1113,err=1113) kk,v1sccor(k,l,i,j) + & ,v2sccor(k,l,i,j) +c write(iout,*) "k=",kk + 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 +C read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j) +C & ,v2sccor(k,l,i,j) 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,*,end=113,err=113) kk,vlor1sccor(k,i,j), + do k=1,nlor_sccor(i,j) + read (isccor,*,end=1113,err=1113) 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(l,i,j)=v0ijsccor + v0sccor(l,-i,j)=v0ijsccor1 + v0sccor(l,i,-j)=v0ijsccor2 + v0sccor(l,-i,-j)=v0ijsccor3 enddo enddo enddo close (isccor) - +#else + write(iout,*) "a tu nie wchodze" + 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=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=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=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,iblock)=v0ijsccor + enddo + enddo + enddo + close (isccor) + +#endif if (lprint) then write (iout,'(/a/)') 'Torsional constants:' do i=1,nsccortyp @@ -646,7 +903,7 @@ C write (iout,*) "Coefficients of the cumulants" endif read (ifourier,*) nloctyp - do i=0,nloctyp-1 + do i=1,nloctyp read (ifourier,*,end=115,err=115) read (ifourier,*,end=115,err=115) (b(ii),ii=1,13) if (lprint) then @@ -655,31 +912,22 @@ C 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)=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) + 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)=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) c CC(1,1,i)=0.0d0 c CC(2,2,i)=0.0d0 c CC(2,1,i)=0.0d0 @@ -688,11 +936,6 @@ c CC(1,2,i)=0.0d0 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) - c Ctilde(1,1,i)=0.0d0 c Ctilde(1,2,i)=0.0d0 c Ctilde(2,1,i)=0.0d0 @@ -701,10 +944,6 @@ c Ctilde(2,2,i)=0.0d0 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) c DD(1,1,i)=0.0d0 c DD(2,2,i)=0.0d0 c DD(2,1,i)=0.0d0 @@ -713,11 +952,6 @@ c DD(1,2,i)=0.0d0 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) - c Dtilde(1,1,i)=0.0d0 c Dtilde(1,2,i)=0.0d0 c Dtilde(2,1,i)=0.0d0 @@ -726,11 +960,6 @@ c Dtilde(2,2,i)=0.0d0 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) - c ee(1,1,i)=1.0d0 c ee(2,2,i)=1.0d0 c ee(2,1,i)=0.0d0 @@ -826,9 +1055,17 @@ C----------------------- LJK potential -------------------------------- endif goto 50 C---------------------- GB or BP potential ----------------------------- - 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), - & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp), - & (alp(i),i=1,ntyp) + 30 do i=1,ntyp + read (isidep,*,end=116,err=116)(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 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), +c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp), +c & (alp(i),i=1,ntyp) C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@ -940,7 +1177,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) C C Define the SC-p interaction constants (hard-coded; old style) C - do i=1,20 + do i=1,ntyp C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates C helix formation) c aad(i,1)=0.3D0*4.0D0**12 @@ -978,7 +1215,7 @@ C if (lprint) then write (iout,*) "Parameters of SC-p interactions:" - do i=1,20 + do i=1,ntyp write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1), & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2) enddo @@ -1027,6 +1264,9 @@ c v3ss=0.0d0 goto 999 113 write (iout,*) "Error reading torsional energy parameters." goto 999 + 1113 write (iout,*) + & "Error reading side-chain torsional energy parameters." + goto 999 114 write (iout,*) "Error reading double torsional energy parameters." goto 999 115 write (iout,*)