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 implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'sizesclu.dat' 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.SCCOR' include 'COMMON.SCROT' character*1 t1,t2,t3 character*1 onelett(4) /"G","A","P","D"/ logical lprint dimension blower(3,3,maxlob) double precision ip,mp C C Body C C Set LPRINT=.TRUE. for debugging dwa16=2.0d0**(1.0d0/6.0d0) lprint=.false. itypro=20 C Assign virtual-bond length vbl=3.8D0 vblinv=1.0D0/vbl vblinv2=vblinv*vblinv #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 #ifdef AIX call flush_(iout) #else call flush(iout) #endif 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 l=1,maxinter do i=1,nsccortyp do j=1,nsccortyp write (iout,*) 'ityp',i,' jtyp',j write (iout,*) 'Fourier constants' do k=1,nterm_sccor(i,j) write (iout,'(2(1pe15.5))') & v1sccor(k,l,i,j),v2sccor(k,l,i,j) enddo write (iout,*) 'Lorenz constants' do k=1,nlor_sccor(i,j) write (iout,'(3(1pe15.5))') & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j) enddo enddo enddo 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.5) 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) 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 50 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 50 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 50 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 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 C C Define the SC-p interaction constants C do i=1,20 do j=1,2 eps_scp(i,j)=-1.5d0 rscp(i,j)=4.0d0 enddo enddo #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 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 C D0CM = 3.78d0 C AKCM = 15.1d0 C AKTH = 11.0d0 C AKCT = 12.0d0 C V1SS =-1.08d0 C V2SS = 7.61d0 C V3SS = 13.7d0 return end