--- /dev/null
+ subroutine parmread(iparm,*)
+C
+C Read the parameters of the probability distributions of the virtual-bond
+C valence angles and the side chains and energy parameters.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.TORSION'
+ include 'COMMON.FFIELD'
+ include 'COMMON.NAMES'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.WEIGHTS'
+ include 'COMMON.ENEPS'
+ include 'COMMON.SCCOR'
+ include 'COMMON.SCROT'
+ include 'COMMON.FREE'
+ character*1 t1,t2,t3
+ character*1 onelett(4) /"G","A","P","D"/
+ logical lprint
+ dimension blower(3,3,maxlob)
+ character*800 controlcard
+ character*256 bondname_t,thetname_t,rotname_t,torname_t,
+ & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
+ & sccorname_t
+ integer ilen
+ external ilen
+ character*16 key
+ integer iparm
+ double precision ip,mp
+C
+C Body
+C
+C Set LPRINT=.TRUE. for debugging
+ dwa16=2.0d0**(1.0d0/6.0d0)
+ lprint=.true.
+ itypro=20
+C Assign virtual-bond length
+ vbl=3.8D0
+ vblinv=1.0D0/vbl
+ vblinv2=vblinv*vblinv
+ call card_concat(controlcard,.true.)
+ wname(4)="WCORRH"
+ do i=1,n_ene
+ key = wname(i)(:ilen(wname(i)))
+ call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
+ enddo
+
+ write (iout,*) "iparm",iparm," myparm",myparm
+c If reading not own parameters, skip assignment
+
+ if (iparm.eq.myparm .or. .not.separate_parset) then
+
+c
+c Setup weights for UNRES
+c
+ wsc=ww(1)
+ wscp=ww(2)
+ welec=ww(3)
+ wcorr=ww(4)
+ wcorr5=ww(5)
+ wcorr6=ww(6)
+ wel_loc=ww(7)
+ wturn3=ww(8)
+ wturn4=ww(9)
+ wturn6=ww(10)
+ wang=ww(11)
+ wscloc=ww(12)
+ wtor=ww(13)
+ wtor_d=ww(14)
+ wvdwpp=ww(16)
+ wstrain=ww(15)
+ wbond=ww(18)
+ wsccor=ww(19)
+
+ endif
+
+ call card_concat(controlcard,.false.)
+
+c Return if not own parameters
+
+ if (iparm.ne.myparm .and. separate_parset) return
+
+ call reads(controlcard,"BONDPAR",bondname_t,bondname)
+ open (ibond,file=bondname_t,status='old')
+ rewind(ibond)
+ call reads(controlcard,"THETPAR",thetname_t,thetname)
+ open (ithep,file=thetname_t,status='old')
+ rewind(ithep)
+ call reads(controlcard,"ROTPAR",rotname_t,rotname)
+ open (irotam,file=rotname_t,status='old')
+ rewind(irotam)
+ call reads(controlcard,"TORPAR",torname_t,torname)
+ open (itorp,file=torname_t,status='old')
+ rewind(itorp)
+ call reads(controlcard,"TORDPAR",tordname_t,tordname)
+ open (itordp,file=tordname_t,status='old')
+ rewind(itordp)
+ call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
+ open (isccor,file=sccorname_t,status='old')
+ rewind(isccor)
+ call reads(controlcard,"FOURIER",fouriername_t,fouriername)
+ open (ifourier,file=fouriername_t,status='old')
+ rewind(ifourier)
+ call reads(controlcard,"ELEPAR",elename_t,elename)
+ open (ielep,file=elename_t,status='old')
+ rewind(ielep)
+ call reads(controlcard,"SIDEPAR",sidename_t,sidename)
+ open (isidep,file=sidename_t,status='old')
+ rewind(isidep)
+ call reads(controlcard,"SCPPAR",scpname_t,scpname)
+ open (iscpp,file=scpname_t,status='old')
+ rewind(iscpp)
+ write (iout,*) "Parameter set:",iparm
+ write (iout,*) "Energy-term weights:"
+ do i=1,n_ene
+ write (iout,'(a16,f10.5)') wname(i),ww(i)
+ enddo
+ write (iout,*) "Sidechain potential file : ",
+ & sidename_t(:ilen(sidename_t))
+#ifndef OLDSCP
+ write (iout,*) "SCp potential file : ",
+ & scpname_t(:ilen(scpname_t))
+#endif
+ write (iout,*) "Electrostatic potential file : ",
+ & elename_t(:ilen(elename_t))
+ write (iout,*) "Cumulant coefficient file : ",
+ & fouriername_t(:ilen(fouriername_t))
+ write (iout,*) "Torsional parameter file : ",
+ & torname_t(:ilen(torname_t))
+ write (iout,*) "Double torsional parameter file : ",
+ & tordname_t(:ilen(tordname_t))
+ write (iout,*) "Backbone-rotamer parameter file : ",
+ & sccorname(:ilen(sccorname))
+ write (iout,*) "Bond & inertia constant file : ",
+ & bondname_t(:ilen(bondname_t))
+ write (iout,*) "Bending parameter file : ",
+ & thetname_t(:ilen(thetname_t))
+ write (iout,*) "Rotamer parameter file : ",
+ & rotname_t(:ilen(rotname_t))
+
+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
+ do i=1,ntyp
+ nbondterm(i)=1
+ read (ibond,*) vbldsc0(1,i),aksc(1,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,*) ijunk,vbldp0,akp,rjunk
+ do i=1,ntyp
+ read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
+ & j=1,nbondterm(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/)')"Force constants virtual bonds:"
+ write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
+ & 'inertia','Pstok'
+ write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
+ do i=1,ntyp
+ write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
+ & vbldsc0(1,i),aksc(1,i),abond0(1,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,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
+ read (ithep,*) (polthet(j,i),j=0,3)
+ read (ithep,*) (gthet(j,i),j=1,3)
+ read (ithep,*) theta0(i),sig0(i),sigc0(i)
+ sigc0(i)=sigc0(i)**2
+ enddo
+ close (ithep)
+ if (lprint) then
+c write (iout,'(a)')
+c & 'Parameters of the virtual-bond valence angles:'
+c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
+c & ' ATHETA0 ',' A1 ',' A2 ',
+c & ' B1 ',' B2 '
+c do i=1,ntyp
+c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
+c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
+c enddo
+c write (iout,'(/a/9x,5a/79(1h-))')
+c & 'Parameters of the expression for sigma(theta_c):',
+c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
+c & ' ALPH3 ',' SIGMA0C '
+c do i=1,ntyp
+c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
+c & (polthet(j,i),j=0,3),sigc0(i)
+c enddo
+c write (iout,'(/a/9x,5a/79(1h-))')
+c & 'Parameters of the second gaussian:',
+c & ' THETA0 ',' SIGMA0 ',' G1 ',
+c & ' G2 ',' G3 '
+c do i=1,ntyp
+c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
+c & sig0(i),(gthet(j,i),j=1,3)
+c enddo
+ 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
+#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,*) nthetyp,ntheterm,ntheterm2,
+ & ntheterm3,nsingle,ndouble
+ nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
+ read (ithep,*) (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)') res1,res2,res3
+ read (ithep,*) aa0thet(i,j,k)
+ read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
+ read (ithep,*)
+ & ((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,*)
+ & (((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
+#endif
+
+#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)') 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,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
+ do j=2,nlob(i)
+ read (irotam,*) bsc(j,i)
+ read (irotam,*) (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
+ write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
+ & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
+c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
+c write (iout,'(a,f10.4,4(16x,f10.4))')
+c & 'Center ',(bsc(j,i),j=1,nlobi)
+c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
+ 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)
+c write (iout,'(a)')
+c do j=1,nlobi
+c ind=0
+c do k=1,3
+c do l=1,k
+c ind=ind+1
+c blower(k,l,j)=gaussc(ind,j,i)
+c enddo
+c enddo
+c enddo
+ 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
+ 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,*)
+ if (i.eq.10) then
+ read (irotam,*)
+ else
+ do j=1,65
+ read(irotam,*) sc_parmin(j,i)
+ enddo
+ endif
+ enddo
+#endif
+ close(irotam)
+#ifdef CRYST_TOR
+C
+C Read torsional parameters in old format
+C
+ read (itorp,*) ntortyp,nterm_old
+ write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
+ read (itorp,*) (itortyp(i),i=1,ntyp)
+ do i=1,ntortyp
+ do j=1,ntortyp
+ read (itorp,'(a)')
+ do k=1,nterm_old
+ read (itorp,*) 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,*) ntortyp
+ read (itorp,*) (itortyp(i),i=1,ntyp)
+ write (iout,*) 'ntortyp',ntortyp
+ do i=1,ntortyp
+ do j=1,ntortyp
+ read (itorp,*) nterm(i,j),nlor(i,j)
+ v0ij=0.0d0
+ si=-1.0d0
+ do k=1,nterm(i,j)
+ read (itorp,*) 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,*) 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)') 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
+ stop "Error in double torsional parameter file"
+ endif
+ read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
+ read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
+ read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
+ read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
+ read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
+ read (itordp,*) ((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,*) nsccortyp
+ read (isccor,*) (isccortyp(i),i=1,ntyp)
+c write (iout,*) 'ntortyp',ntortyp
+ maxinter=3
+cc maxinter is maximum interaction sites
+ do l=1,maxinter
+ do i=1,nsccortyp
+ do j=1,nsccortyp
+ read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
+ v0ijsccor=0.0d0
+ si=-1.0d0
+
+ do k=1,nterm_sccor(i,j)
+ read (isccor,*) kk,v1sccor(k,l,i,j)
+ & ,v2sccor(k,l,i,j)
+ v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
+ si=-si
+ enddo
+ do k=1,nlor_sccor(i,j)
+ read (isccor,*) kk,vlor1sccor(k,i,j),
+ & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+ v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
+ &(1+vlor3sccor(k,i,j)**2)
+ enddo
+ v0sccor(i,j)=v0ijsccor
+ enddo
+ enddo
+ enddo
+ close (isccor)
+
+ if (lprint) then
+ write (iout,'(/a/)') 'Torsional constants:'
+ do i=1,nsccortyp
+ do j=1,nsccortyp
+ write (iout,*) 'ityp',i,' jtyp',j
+ write (iout,*) 'Fourier constants'
+ do k=1,nterm_sccor(i,j)
+ write(iout,'(2(1pe15.5))')(v1sccor(k,l,i,j),v2sccor(k,l,i,j),
+ & l=1,maxinter)
+ enddo
+ write (iout,*) 'Lorenz constants'
+ do k=1,nlor_sccor(i,j)
+ write (iout,'(3(1pe15.5))')
+ & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+ enddo
+ enddo
+ enddo
+ endif
+
+C
+C 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,*) nloctyp
+ do i=1,nloctyp
+ read (ifourier,*)
+ read (ifourier,*) (b(ii,i),ii=1,13)
+ if (lprint) then
+ write (iout,*) 'Type',i
+ write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
+ endif
+ 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)
+ 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)
+ 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)
+ enddo
+ if (lprint) then
+ do i=1,nloctyp
+ write (iout,*) 'Type',i
+ write (iout,*) 'B1'
+c write (iout,'(f10.5)') B1(:,i)
+ write(iout,*) B1(1,i),B1(2,i)
+ write (iout,*) 'B2'
+c write (iout,'(f10.5)') B2(:,i)
+ 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,'(/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,*) ((epp(i,j),j=1,2),i=1,2)
+ read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
+ read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
+ read (ielep,*) ((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,*) ipot,expon
+ if (ipot.lt.1 .or. ipot.gt.6) then
+ write (iout,'(2a)') 'Error while reading SC interaction',
+ & 'potential file - unknown potential type.'
+ stop
+ endif
+ expon2=expon/2
+ 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,*)((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,*)((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,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
+ & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(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)=(chip0(i)-1.0D0)/(chip0(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,*)((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
+
+ 60 continue
+ close (isidep)
+C-----------------------------------------------------------------------
+C Calculate the "working" parameters of SC interactions.
+
+ 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
+ write (iout,*) "IPOT=",ipot
+ 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
+
+C
+C Define the SC-p interaction constants
+C
+#ifdef OLDSCP
+ 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,*) (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
+
+ if (lprint) 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
+ end
projects / unres.git / blobdiff