X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_Eshel%2Fparmread.F;fp=source%2Funres%2Fsrc_Eshel%2Fparmread.F;h=cd9bda2c8b24d09b74b97c4effd8558e21e31e25;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_Eshel/parmread.F b/source/unres/src_Eshel/parmread.F new file mode 100644 index 0000000..cd9bda2 --- /dev/null +++ b/source/unres/src_Eshel/parmread.F @@ -0,0 +1,1036 @@ + 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"/ + 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),j=1,2), + & (bthet(j,i),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 + 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),j=1,2),(bthet(j,i),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),j=1,2),(10*bthet(j,i),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),j=1,2), + & (bthet(j,i),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 + write (2,*) "End reading THETA_PDB" + 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) + 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) + 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) + 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) +c write (iout,*) 'ntortyp',ntortyp + do i=1,ntortyp + do j=1,ntortyp + read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j) + v0ij=0.0d0 + si=-1.0d0 + do k=1,nterm(i,j) + read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j) + v0ij=v0ij+si*v1(k,i,j) + si=-si + enddo + do k=1,nlor(i,j) + 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)=v0ij + 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) + write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j) + enddo + write (iout,*) 'Lorenz constants' + do k=1,nlor(i,j) + 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 i=1,ntortyp + do j=1,ntortyp + do k=1,ntortyp + 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 + 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), + & 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 + if (lprint) then + write (iout,*) + write (iout,*) 'Constants for double torsionals' + do i=1,ntortyp + do j=1,ntortyp + do k=1,ntortyp + write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k, + & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k) + 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) + 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)) + 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)) + enddo + write (iout,*) + 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=1113,err=1113) nsccortyp + read (isccor,*,end=1113,err=1113) (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=1113,err=1113) + & nterm_sccor(i,j),nlor_sccor(i,j) + v0ijsccor=0.0d0 + si=-1.0d0 + + do k=1,nterm_sccor(i,j) + read (isccor,*,end=1113,err=1113) 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=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 + enddo + enddo + enddo + close (isccor) + + if (lprint) then + write (iout,'(/a/)') 'Torsional constants:' + do i=1,nsccortyp + do j=1,nsccortyp + write (iout,*) 'ityp',i,' jtyp',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 + endif +C +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=1,nloctyp + 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) +c b1(1,i)=0.0d0 +c b1(2,i)=0.0d0 + 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) +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) +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) +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) +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) +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) +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 + 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,*) + & "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