+++ /dev/null
- 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
projects / unres.git / blobdiff