X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_CSA_DiL%2Fparmread.F;fp=source%2Funres%2Fsrc_CSA_DiL%2Fparmread.F;h=0000000000000000000000000000000000000000;hp=44d0370eb440288bf544e2dc65934adedb59e559;hb=2a226bfc86eabc6e4eae0c3ad1cbc3cb5417a05a;hpb=a0e685f844163003749ba91dfbf4644bcc8cfa30 diff --git a/source/unres/src_CSA_DiL/parmread.F b/source/unres/src_CSA_DiL/parmread.F deleted file mode 100644 index 44d0370..0000000 --- a/source/unres/src_CSA_DiL/parmread.F +++ /dev/null @@ -1,1132 +0,0 @@ - subroutine parmread -C -C Read the parameters of the probability distributions of the virtual-bond -C valence angles and the side chains and energy parameters. -C -C Important! Energy-term weights ARE NOT read here; they are read from the -C main input file instead, because NO defaults have yet been set for these -C parameters. -C - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' -#ifdef MPI - include "mpif.h" - integer IERROR -#endif - include 'COMMON.IOUNITS' - include 'COMMON.CHAIN' - include 'COMMON.INTERACT' - include 'COMMON.GEO' - include 'COMMON.LOCAL' - include 'COMMON.TORSION' - include 'COMMON.SCCOR' - include 'COMMON.SCROT' - include 'COMMON.FFIELD' - include 'COMMON.NAMES' - include 'COMMON.SBRIDGE' - include 'COMMON.MD_' - 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) - character*3 lancuch,ucase -C -C For printing parameters after they are read set the following in the UNRES -C C-shell script: -C -C setenv PRINT_PARM YES -C -C To print parameters in LaTeX format rather than as ASCII tables: -C -C setenv LATEX YES -C - call getenv_loc("PRINT_PARM",lancuch) - lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y") - call getenv_loc("LATEX",lancuch) - LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y") -C - dwa16=2.0d0**(1.0d0/6.0d0) - itypro=20 -C Assign virtual-bond length - vbl=3.8D0 - vblinv=1.0D0/vbl - vblinv2=vblinv*vblinv -c -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 - do i=1,ntyp - nbondterm(i)=1 - read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i) - dsc(i) = vbldsc0(1,i) - if (i.eq.10) then - dsc_inv(i)=0.0D0 - else - dsc_inv(i)=1.0D0/dsc(i) - endif - enddo -#else - read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok - do i=1,ntyp - 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) - if (i.eq.10) then - dsc_inv(i)=0.0D0 - else - dsc_inv(i)=1.0D0/dsc(i) - endif - enddo -#endif - if (lprint) then - write(iout,'(/a/)')"Dynamic constants of the interaction sites:" - write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass', - & 'inertia','Pstok' - write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok - do i=1,ntyp - write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i), - & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i) - do j=2,nbondterm(i) - write (iout,'(13x,3f10.5)') - & vbldsc0(j,i),aksc(j,i),abond0(j,i) - enddo - enddo - endif -#ifdef CRYST_THETA -C -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,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 - 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 - write (iout,'(a)') - & 'Parameters of the virtual-bond valence angles:' - write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:', - & ' ATHETA0 ',' A1 ',' A2 ', - & ' B1 ',' B2 ' - do i=1,ntyp - write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i, - & 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):', - & ' ALPH0 ',' ALPH1 ',' ALPH2 ', - & ' ALPH3 ',' SIGMA0C ' - do i=1,ntyp - write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i, - & (polthet(j,i),j=0,3),sigc0(i) - enddo - write (iout,'(/a/9x,5a/79(1h-))') - & 'Parameters of the second gaussian:', - & ' THETA0 ',' SIGMA0 ',' G1 ', - & ' G2 ',' G3 ' - do i=1,ntyp - write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i), - & sig0(i),(gthet(j,i),j=1,3) - enddo - else - write (iout,'(a)') - & 'Parameters of the virtual-bond valence angles:' - write (iout,'(/a/9x,5a/79(1h-))') - & 'Coefficients of expansion', - & ' theta0 ',' a1*10^2 ',' a2*10^2 ', - & ' 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,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):', - & ' alpha0 ',' alph1 ',' alph2 ', - & ' alhp3 ',' sigma0c ' - do i=1,ntyp - write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i), - & (polthet(j,i),j=0,3),sigc0(i) - enddo - write (iout,'(/a/9x,5a/79(1h-))') - & 'Parameters of the second gaussian:', - & ' theta0 ',' sigma0*10^2 ',' G1*10^-1', - & ' G2 ',' G3*10^1 ' - do i=1,ntyp - write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i), - & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0 - enddo - endif - endif -#else -C -C Read the parameters of Utheta determined from ab initio surfaces -C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 -C - read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2, - & 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 l=1,ntheterm - aathet(l,i,j,k)=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 - 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 - 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) - 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) - 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), - & 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 - 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) - enddo - aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1) - aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j) - enddo - do l=1,ntheterm - aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1) - enddo - aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1) - enddo -C -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 - 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)') - & (l,aathet(l,i,j,k),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) - enddo - enddo - do l=1,ntheterm3 - write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') - & "f+",l,"f-",l,"g+",l,"g-",l - 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) - enddo - enddo - enddo - enddo - enddo - enddo - call flush(iout) - endif - write (2,*) "Start reading THETA_PDB" - do i=1,ntyp - 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 - write (2,*) "End reading THETA_PDB" - enddo - close (ithep_pdb) -#endif - close(ithep) -#ifdef CRYST_SC -C -C Read the parameters of the probability distribution/energy expression -C of the side chains. -C - do i=1,ntyp - read (irotam,'(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,*,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 - 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 - 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 - endif - enddo - close (irotam) - if (lprint) then - write (iout,'(/a)') 'Parameters of side-chain local geometry' - do i=1,ntyp - nlobi=nlob(i) - if (nlobi.gt.0) then - if (LaTeX) then - write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i), - & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i) - write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') - & 'log h',(bsc(j,i),j=1,nlobi) - write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') - & 'x',((censc(k,j,i),k=1,3),j=1,nlobi) - do k=1,3 - write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') - & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi) - enddo - else - write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi) - write (iout,'(a,f10.4,4(16x,f10.4))') - & 'Center ',(bsc(j,i),j=1,nlobi) - write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3), - & j=1,nlobi) - write (iout,'(a)') - endif - endif - enddo - endif -#else -C -C Read scrot parameters for potentials determined from all-atom AM1 calculations -C added by Urszula Kozlowska 07/11/2007 -C - do i=1,ntyp - read (irotam,*,end=112,err=112) - if (i.eq.10) then - read (irotam,*,end=112,err=112) - else - do j=1,65 - read(irotam,*,end=112,err=112) sc_parmin(j,i) - enddo - endif - enddo -C -C Read the parameters of the probability distribution/energy expression -C of the side chains. -C - 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) -#endif - close(irotam) - -#ifdef CRYST_TOR -C -C Read torsional parameters in old format -C - read (itorp,*,end=113,err=113) ntortyp,nterm_old - if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old - read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp) - do i=1,ntortyp - do j=1,ntortyp - read (itorp,'(a)') - do k=1,nterm_old - read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i) - enddo - enddo - enddo - close (itorp) - if (lprint) then - write (iout,'(/a/)') 'Torsional constants:' - do i=1,ntortyp - do j=1,ntortyp - write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old) - write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old) - enddo - enddo - endif -#else -C -C Read torsional parameters -C - read (itorp,*,end=113,err=113) ntortyp - read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp) - do iblock=1,2 - do i=-ntyp,-1 - itortyp(i)=-itortyp(-i) - enddo -c write (iout,*) 'ntortyp',ntortyp - do i=0,ntortyp-1 - 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) - 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,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) - read (itorp,*,end=113,err=113) kk,vlor1(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 - v0(-i,-j,iblock)=v0ij - enddo - enddo - enddo - close (itorp) - if (lprint) then - 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,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,iblock) - write (iout,'(3(1pe15.5))') - & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) - enddo - enddo - enddo - endif -C -C 6/23/01 Read parameters for double torsionals -C - 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 -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 - call MPI_Finalize(Ierror) -#endif - stop "Error in double torsional parameter file" - endif - 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 Matrix 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 Matrix 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,*) 'Constants for double torsionals' - do iblock=1,2 - do i=1,ntortyp - do j=-ntortyp,ntortyp - do k=-ntortyp,ntortyp - write (iout,*) 'ityp',i,' jtyp',j,' ktyp',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,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) - enddo - write (iout,*) - write (iout,*) 'Pairs of angles:' - 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,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 -C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local -C correlation energies. -C - read (isccor,*,end=119,err=119) nterm_sccor - do i=1,20 - do j=1,20 - read (isccor,'(a)') - do k=1,nterm_sccor - read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j), - & v2sccor(k,i,j) - enddo - enddo - enddo - close (isccor) - if (lprint) then - write (iout,'(/a/)') 'Torsional constants of SCCORR:' - do i=1,20 - do j=1,20 - write (iout,*) 'ityp',i,' jtyp',j - do k=1,nterm_sccor - write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j) - enddo - enddo - enddo - endif -C -C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local -C interaction energy of the Gly, Ala, and Pro prototypes. -C - if (lprint) then - write (iout,*) - 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) - if (lprint) then - write (iout,*) 'Type',i - write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),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)=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)=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 -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) -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) -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) -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) -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 - if (lprint) then - do i=1,nloctyp - write (iout,*) 'Type',i - write (iout,*) 'B1' - write(iout,*) B1(1,i),B1(2,i) - write (iout,*) 'B2' - write(iout,*) B2(1,i),B2(2,i) - write (iout,*) 'CC' - do j=1,2 - write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i) - enddo - write(iout,*) 'DD' - do j=1,2 - write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i) - enddo - write(iout,*) 'EE' - do j=1,2 - write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i) - enddo - enddo - endif -C -C Read electrostatic-interaction parameters -C - if (lprint) then - write (iout,*) - write (iout,'(/a)') 'Electrostatic interaction constants:' - write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') - & 'IT','JT','APP','BPP','AEL6','AEL3' - endif - read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2) - read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2) - read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2) - read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2) - close (ielep) - do i=1,2 - do j=1,2 - rri=rpp(i,j)**6 - app (i,j)=epp(i,j)*rri*rri - 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 - if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j), - & ael6(i,j),ael3(i,j) - enddo - enddo -C -C Read side-chain interaction parameters. -C - read (isidep,*,end=117,err=117) ipot,expon - if (ipot.lt.1 .or. ipot.gt.5) then - write (iout,'(2a)') 'Error while reading SC interaction', - & 'potential file - unknown potential type.' -#ifdef MPI - call MPI_Finalize(Ierror) -#endif - stop - endif - expon2=expon/2 - if(me.eq.king) - & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot), - & ', 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), - & (sigma0(i),i=1,ntyp) - if (lprint) then - write (iout,'(/a/)') 'Parameters of the LJ potential:' - write (iout,'(a/)') 'The epsilon array:' - call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) - write (iout,'(/a)') 'One-body parameters:' - write (iout,'(a,4x,a)') 'residue','sigma' - write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp) - endif - goto 50 -C----------------------- LJK potential -------------------------------- - 20 read (isidep,*,end=116,err=116)((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:' - write (iout,'(a/)') 'The epsilon array:' - call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) - write (iout,'(/a)') 'One-body parameters:' - write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 ' - write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i), - & i=1,ntyp) - 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) -C For the GB potential convert sigma'**2 into chi' - if (ipot.eq.4) then - do i=1,ntyp - chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0) - enddo - endif - if (lprint) then - write (iout,'(/a/)') 'Parameters of the BP potential:' - write (iout,'(a/)') 'The epsilon array:' - call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) - write (iout,'(/a)') 'One-body parameters:' - write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2', - & ' chip ',' alph ' - write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i), - & chip(i),alp(i),i=1,ntyp) - endif - goto 50 -C--------------------- GBV potential ----------------------------------- - 40 read (isidep,*,end=116,err=116)((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 - write (iout,'(/a/)') 'Parameters of the GBV potential:' - write (iout,'(a/)') 'The epsilon array:' - call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) - write (iout,'(/a)') 'One-body parameters:' - write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ', - & 's||/s_|_^2',' chip ',' alph ' - write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i), - & sigii(i),chip(i),alp(i),i=1,ntyp) - endif - 50 continue - close (isidep) -C----------------------------------------------------------------------- -C Calculate the "working" parameters of SC interactions. - do i=2,ntyp - do j=1,i-1 - eps(i,j)=eps(j,i) - enddo - enddo - do i=1,ntyp - do j=i,ntyp - sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2) - sigma(j,i)=sigma(i,j) - rs0(i,j)=dwa16*sigma(i,j) - rs0(j,i)=rs0(i,j) - enddo - enddo - if (lprint) write (iout,'(/a/10x,7a/72(1h-))') - & 'Working parameters of the SC interactions:', - & ' a ',' b ',' augm ',' sigma ',' r0 ', - & ' chi1 ',' chi2 ' - do i=1,ntyp - do j=i,ntyp - epsij=eps(i,j) - if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then - rrij=sigma(i,j) - else - rrij=rr0(i)+rr0(j) - endif - r0(i,j)=rrij - r0(j,i)=rrij - rrij=rrij**expon - 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) - if (ipot.gt.2) then - sigt1sq=sigma0(i)**2 - sigt2sq=sigma0(j)**2 - sigii1=sigii(i) - sigii2=sigii(j) - ratsig1=sigt2sq/sigt1sq - ratsig2=1.0D0/ratsig1 - chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1) - if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2) - rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq) - else - rsum_max=sigma(i,j) - endif -c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then - sigmaii(i,j)=rsum_max - sigmaii(j,i)=rsum_max -c else -c sigmaii(i,j)=r0(i,j) -c sigmaii(j,i)=r0(i,j) -c endif -cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max - if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then - r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij - augm(i,j)=epsij*r_augm**(2*expon) -c augm(i,j)=0.5D0**(2*expon)*aa(i,j) - augm(j,i)=augm(i,j) - else - augm(i,j)=0.0D0 - augm(j,i)=0.0D0 - 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), - & sigma(i,j),r0(i,j),chi(i,j),chi(j,i) - endif - enddo - enddo -#ifdef OLDSCP -C -C Define the SC-p interaction constants (hard-coded; old style) -C - do i=1,20 -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 -C Following line for constants currently implemented -C "Hard" SC-p repulsion (gives correct turn spacing in helices) - aad(i,1)=1.5D0*4.0D0**12 -c aad(i,1)=0.17D0*5.6D0**12 - aad(i,2)=aad(i,1) -C "Soft" SC-p repulsion - bad(i,1)=0.0D0 -C Following line for constants currently implemented -c aad(i,1)=0.3D0*4.0D0**6 -C "Hard" SC-p repulsion - bad(i,1)=3.0D0*4.0D0**6 -c bad(i,1)=-2.0D0*0.17D0*5.6D0**6 - bad(i,2)=bad(i,1) -c aad(i,1)=0.0D0 -c aad(i,2)=0.0D0 -c bad(i,1)=1228.8D0 -c bad(i,2)=1228.8D0 - enddo -#else -C -C 8/9/01 Read the SC-p interaction constants from file -C - do i=1,ntyp - read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2) - enddo - do i=1,ntyp - aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12 - aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12 - 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 - - if (lprint) then - write (iout,*) "Parameters of SC-p interactions:" - do i=1,20 - 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 - endif -#endif -C -C Define the constants of the disulfide bridge -C - ebr=-5.50D0 -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 -c akcm=0.0d0 -c akth=0.0d0 -c akct=0.0d0 -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 - return - 111 write (iout,*) "Error reading bending energy parameters." - goto 999 - 112 write (iout,*) "Error reading rotamer energy parameters." - goto 999 - 113 write (iout,*) "Error reading torsional energy parameters." - goto 999 - 114 write (iout,*) "Error reading double torsional energy parameters." - goto 999 - 115 write (iout,*) - & "Error reading cumulant (multibody energy) parameters." - goto 999 - 116 write (iout,*) "Error reading electrostatic energy parameters." - goto 999 - 117 write (iout,*) "Error reading side chain interaction parameters." - goto 999 - 118 write (iout,*) "Error reading SCp interaction parameters." - goto 999 - 119 write (iout,*) "Error reading SCCOR parameters" - 999 continue -#ifdef MPI - call MPI_Finalize(Ierror) -#endif - stop - return - end - - - subroutine getenv_loc(var, val) - character(*) var, val - -#ifdef WINIFL - character(2000) line - external ilen - - open (196,file='env',status='old',readonly,shared) - iread=0 -c write(*,*)'looking for ',var -10 read(196,*,err=11,end=11)line - iread=index(line,var) -c write(*,*)iread,' ',var,' ',line - if (iread.eq.0) go to 10 -c write(*,*)'---> ',line -11 continue - if(iread.eq.0) then -c write(*,*)'CHUJ' - val='' - else - iread=iread+ilen(var)+1 - read (line(iread:),*,err=12,end=12) val -c write(*,*)'OK: ',var,' = ',val - endif - close(196) - return -12 val='' - close(196) -#elif (defined CRAY) - integer lennam,lenval,ierror -c -c getenv using a POSIX call, useful on the T3D -c Sept 1996, comment out error check on advice of H. Pritchard -c - lennam = len(var) - if(lennam.le.0) stop '--error calling getenv--' - call pxfgetenv(var,lennam,val,lenval,ierror) -c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--' -#else - call getenv(var,val) -#endif - - return - end