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) 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 END DO 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) END DO v0sccor(i,j) = v0ijsccor END DO END DO END DO 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) END DO 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) END DO END DO END DO END IF 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" END IF 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) END IF 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.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 GB or 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),sigmap1(i,j),sigmap2(i,j), & chis(i,j),chis(j,i), & nstate(i,j),(wstate(k,i,j),k=1,4), & dhead(1,1,i,j), & dhead(1,2,i,j), & dhead(2,1,i,j), & dhead(2,2,i,j), & dtail(1,i,j),dtail(2,i,j), & 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),sigiso1(i,j),sigiso2(i,j),epsintab(i,j) END DO END DO c! write (*,*) "nstate gly-gly", nstate(10,10) 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 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) sigmap1(i,j) = sigmap1(j,i) sigmap2(i,j) = sigmap2(j,i) sigiso1(i,j) = sigiso1(j,i) sigiso2(i,j) = sigiso2(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 dhead(2,1,i,j) = dhead(1,1,j,i) dhead(2,2,i,j) = dhead(1,2,j,i) dhead(1,1,i,j) = dhead(2,1,j,i) dhead(1,2,i,j) = dhead(2,2,j,i) dtail(1,i,j) = dtail(1,j,i) dtail(2,i,j) = dtail(2,j,i) c! DO k = 1, 2 c! DO l = 1, 2 c! write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i) c! write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j) c! dhead(l,k,i,j) = dhead(k,l,j,i) c! END DO c! 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 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), & sigmap1(i,j),sigmap2(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(1,i,j),dtail(2,i,j), & 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),sigiso1(i,j),sigiso2(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) 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).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) 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 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),sigmap1(i,j),sigmap2(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(1,i,j),dtail(2,i,j), & 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),sigiso1(i,j),sigiso2(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,*) "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