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 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 read (ifourier,*,end=115,err=115) nloctyp do i=1,nloctyp read (ifourier,*) read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13) #ifdef NEWCORR read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3) read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3) read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1) read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1) read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1) #endif if (lprint) then write (iout,*) 'Type',i write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13) do ii=1,3 write (iout,'("Bnew1",i3,2f10.5)') ii,(bnew1(ii,k,i),k=1,2) write (iout,'("Bnew2",i3,2f10.5)') ii,(bnew2(ii,k,i),k=1,2) enddo do ii=1,1 write (iout,'("EEnew",i3,2f10.5)') ii,eenew(ii,i) enddo endif #ifndef NEWCORR B1(1,i) = b(3,i) B1(2,i) = b(5,i) B1tilde(1,i) = b(3,i) B1tilde(2,i) =-b(5,i) B2(1,i) = b(2,i) B2(2,i) = b(4,i) #endif CC(1,1,i)= b(7,i) CC(2,2,i)=-b(7,i) CC(2,1,i)= b(9,i) CC(1,2,i)= b(9,i) Ctilde(1,1,i)=b(7,i) Ctilde(1,2,i)=b(9,i) Ctilde(2,1,i)=-b(9,i) Ctilde(2,2,i)=b(7,i) DD(1,1,i)= b(6,i) DD(2,2,i)=-b(6,i) DD(2,1,i)= b(8,i) DD(1,2,i)= b(8,i) Dtilde(1,1,i)=b(6,i) Dtilde(1,2,i)=b(8,i) Dtilde(2,1,i)=-b(8,i) Dtilde(2,2,i)=b(6,i) #ifdef NEWCORR EEold(1,1,i)= b(10,i)+b(11,i) EEold(2,2,i)=-b(10,i)+b(11,i) EEold(2,1,i)= b(12,i)-b(13,i) EEold(1,2,i)= b(12,i)+b(13,i) c EEold(1,1,-i)= b(10,i)+b(11,i) c EEold(2,2,-i)=-b(10,i)+b(11,i) c EEold(2,1,-i)=-b(12,i)+b(13,i) c EEold(1,2,-i)=-b(12,i)-b(13,i) #else EE(1,1,i)= b(10,i)+b(11,i) EE(2,2,i)=-b(10,i)+b(11,i) EE(2,1,i)= b(12,i)-b(13,i) EE(1,2,i)= b(12,i)+b(13,i) #endif 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.6) 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,50) 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 60 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 60 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 60 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 goto 60 C--------------------- Momo potential ----------------------------------- 50 continue read (isidep,*) (icharge(i),i=1,ntyp) c write (2,*) "icharge",(icharge(i),i=1,ntyp) do i=1,ntyp do j=1,i c! write (*,*) "Im in ", i, " ", j read(isidep,*) & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & (alphasur(k,i,j),k=1,4),sigmap(i,j),sigmap(j,i), & chis(i,j),chis(j,i), & nstate(i,j),(wstate(k,i,j),k=1,4), & ((dhead(l,k,i,j),k=1,2),l=1,2),dtail(i,j),dtail(j,i), & epshead(i,j),sig0head(i,j), & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & alphapol(i,j),alphapol(j,i), & (alphiso(k,i,j),k=1,4),sigiso(i,j),sigiso(j,i),epsintab(i,j) END DO END DO c! write (*,*) "Parameters read in" c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TRY TO TAKE PARAMETERS c! FROM HIS-ARG. c! c! DISABLE IT AT >>YOUR OWN PERIL<< c! DO i = 1, ntyp DO j = i+1, ntyp eps(i,j) = eps(j,i) sigma(i,j) = sigma(j,i) nstate(i,j) = nstate(j,i) DO k = 1, 4 alphasur(k,i,j) = alphasur(k,j,i) wstate(k,i,j) = wstate(k,j,i) alphiso(k,i,j) = alphiso(k,j,i) END DO DO k = 1, 2 DO l = 1, 2 dhead(k,l,i,j) = dhead(k,l,j,i) END DO END DO epshead(i,j) = epshead(j,i) sig0head(i,j) = sig0head(j,i) DO k = 1, 2 wqdip(k,i,j) = wqdip(k,j,i) END DO wquad(i,j) = wquad(j,i) epsintab(i,j) = epsintab(j,i) END DO END DO c! IF (1.eq.0) THEN c! write (*,*) "" c! write (*,*) "" c! write (*,'(f17.15)') ((eps(i,j),j=i,ntyp), i=1, ntyp) c! write (*,'(f17.15)') (sigma0(i),i=1,ntyp) c! write (*,'(f17.15)') (sigii(i),i=1,ntyp) c! write (*,'(f17.15)') (chip(i),i=1,ntyp) c! write (*,'(f17.15)') (alp(i),i=1,ntyp) c! write (*,*) "" c! write (*,*) "" c! END IF c! This is older parameter-filling loop c! DO i = 1, ntyp c! DO j = 1, ntyp c c IF ((sig0head(i,j).EQ.0.0d0).AND.(sig0head(j,i).NE.0.0d0)) THEN c sig0head(i,j) = sig0head(j,i) c ELSE IF ((sig0head(j,i).EQ.0.0d0).AND.(sig0head(i,j).NE.0.0d0)) c & THEN c sig0head(j,i) = sig0head(i,j) c END IF c c IF ((epshead(i,j).EQ.0.0d0).AND.(epshead(j,i).NE.0.0d0)) THEN c epshead(i,j) = epshead(j,i) c ELSE IF ((epshead(j,i).EQ.0.0d0).AND.(epshead(i,j).NE.0.0d0)) c & THEN c epshead(j,i) = epshead(i,j) c END IF c c IF ((wquad(i,j).EQ.0.0d0).AND.(wquad(j,i).NE.0.0d0)) THEN c wquad(i,j) = wquad(j,i) c ELSE IF ((wquad(j,i).EQ.0.0d0).AND.(wquad(i,j).NE.0.0d0)) c & THEN c wquad(j,i) = wquad(i,j) c END IF c c IF ((sigiso(i,j).EQ.0.0d0).AND.sigiso(j,i).NE.0.0d0)) THEN c sigiiso(i,j) = sigiso(j,i) c ELSE IF ((epshead(j,i).EQ.0.0d0).AND.(epshead(i,j).NE.0.0d0)) c & THEN c sigiso(j,i) = sigiso(i,j) c END IF c! DO k = 1, 4 c! IF ((wstate(k,i,j).EQ.0.0d0).AND. c! & (wstate(k,j,i).NE.0.0d0)) THEN c! wstate(k,i,j) = wstate(k,j,i) c! ELSE IF ((wstate(k,j,i).EQ.0.0d0).AND. c! & (wstate(k,i,j).NE.0.0d0)) THEN c! wstate(k,j,i) = wstate(k,i,j) c! END IF c! END DO c! END DO c! END DO c! THE LOOP MAKING MOMO_TABLES SYMMETRIC ENDS HERE if (.not.lprint) goto 70 write (iout,'(a)') & "Parameters of the new physics-based SC-SC interaction potential" write (iout,'(/7a)') 'Residues',' epsGB',' rGB', & ' chi1GB',' chi2GB',' chip1GB',' chip2GB' do i=1,ntyp do j=1,i write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))') & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i), & chipp(i,j),chipp(j,i) enddo enddo write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2', & ' alphasur3',' alphasur4',' sigmap1',' sigmap2', & ' chis1',' chis2' do i=1,ntyp do j=1,i write (iout,'(2(a3,1x),8(0pf10.3))') & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4), & sigmap(i,j),sigmap(j,i),chis(i,j),chis(j,i) enddo enddo write (iout,'(/14a)') 'Residues',' nst ',' wst1', & ' wst2',' wst3',' wst4',' dh11',' dh21', & ' dh12',' dh22',' dt1',' dt2',' epsh1', & ' sigh' do i=1,ntyp do j=1,i write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)') & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4), & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(i,j),dtail(j,i), & epshead(i,j),sig0head(i,j) enddo enddo write (iout,'(/12a)') 'Residues',' ch1',' ch2', & ' rborn1',' rborn2',' wqdip1',' wqdip2', & ' wquad' do i=1,ntyp do j=1,i write (iout,'(2(a3,1x),2i4,5f10.3)') & restyp(i),restyp(j),icharge(i),icharge(j), & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j) enddo enddo write (iout,'(/12a)') 'Residues', & ' alphpol1', & ' alphpol2',' alphiso1',' alpiso2', & ' alpiso3',' alpiso4',' sigiso1',' sigiso2', & ' epsin' do i=1,ntyp do j=1,i write (iout,'(2(a3,1x),11f10.3)') & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i), & (alphiso(k,i,j),k=1,4),sigiso(i,j),sigiso(j,i),epsintab(i,j) enddo enddo goto 70 C For the GB potential convert sigma'**2 into chi' DO i=1,ntyp chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0) END DO 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) END IF GOTO 60 60 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) do k=1,4 alphasur(k,j,i)=alphasur(k,i,j) alphiso(k,j,i)=alphiso(k,i,j) wstate(k,j,i)=wstate(k,i,j) enddo do k=1,2 wqdip(k,j,i)=wqdip(k,i,j) enddo do k=1,2 do l=1,2 dhead(l,k,j,i)=dhead(l,k,i,j) enddo enddo enddo enddo IF (ipot.LT.6) THEN 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 END IF 70 continue 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 .or. ipot.eq.6) & THEN rrij=sigma(i,j) ELSE rrij=rr0(i)+rr0(j) END IF 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).AND.(ipot.LT.6)) 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) END IF 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 if (ipot.lt.6) 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) else write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4, & i3,40f10.4)') & 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), & icharge(i),icharge(j),chipp(i,j),chipp(j,i), & (alphasur(k,i,j),k=1,4),sigmap(i,j),sigmap(j,i), & chis(i,j),chis(j,i), & nstate(i,j),(wstate(k,i,j),k=1,4), & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(i,j),dtail(j,i), & epshead(i,j),sig0head(i,j), & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j), & alphapol(i,j),alphapol(j,i), & (alphiso(k,i,j),k=1,4),sigiso(i,j) endif 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,'(/a)') "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 c! print *, " " c! print *, "END OF PARMREAD" c! print *, "eps = ", eps(9,9), "sigma = ", sigma(9,9) c! print *, "chi1 = ", chi(9,9), "chi2 = ", chi(9,9) c! print *, "chip1 = ", chipp(9,9), "chip2 = ", chipp(9,9) c! print *, "sig1 = ", sigmap(9,9), "sig2 = ", sigmap(9,9) c! print *, "chis1 = ", chis(9,9)," chis2 = ", chis(9,9) c! print *, "END OF PARMREAD" c! print *, " " 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