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 none 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' #ifdef FIVEDIAG include 'COMMON.LAGRANGE.5diag' #else include 'COMMON.LAGRANGE' #endif #ifdef LANG0 #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #else include 'COMMON.LANGEVIN' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.SHIELD' 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 double precision blower(3,3,maxlob) character*3 string character*3 lancuch,ucase character*1000 weightcard character*4 res1 integer i,ii,j,jj,k,kk,l,ll,lll,llll,m,mm,n,iblock,junk,ijunk, & nkcctyp,maxinter double precision akl,v0ij,si,rri,epsij,v0ijsccor,epsijlip,rjunk, & sigt2sq,sigt1sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm, & rrij,sigeps double precision dwa16 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 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,vbldpdum,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,vbldpdum,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 C reading lipid parameters if (lprint) then write (iout,*) "iliptranpar",iliptranpar call flush(iout) endif read(iliptranpar,*) pepliptran do i=1,ntyp read(iliptranpar,*) liptranene(i) enddo close(iliptranpar) #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 IF (tor_mode.eq.0) THEN 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=-ntyp1,-1 ithetyp(i)=-ithetyp(-i) enddo do iblock=1,2 do i=-maxthetyp,maxthetyp do j=-maxthetyp,maxthetyp do k=-maxthetyp,maxthetyp aa0thet(i,j,k,iblock)=0.0d0 do l=1,ntheterm aathet(l,i,j,k,iblock)=0.0d0 enddo do l=1,ntheterm2 do m=1,nsingle bbthet(m,l,i,j,k,iblock)=0.0d0 ccthet(m,l,i,j,k,iblock)=0.0d0 ddthet(m,l,i,j,k,iblock)=0.0d0 eethet(m,l,i,j,k,iblock)=0.0d0 enddo enddo do l=1,ntheterm3 do m=1,ndouble do mm=1,ndouble ffthet(mm,m,l,i,j,k,iblock)=0.0d0 ggthet(mm,m,l,i,j,k,iblock)=0.0d0 enddo enddo enddo enddo enddo enddo enddo c VAR:iblock means terminally blocking group 1=non-proline 2=proline do iblock=1,2 c VAR:ntethtyp is type of theta potentials type currently 0=glycine c VAR:1=non-glicyne non-proline 2=proline c VAR:negative values for D-aminoacid do i=0,nthetyp do j=-nthetyp,nthetyp do k=-nthetyp,nthetyp read (ithep,'(6a)',end=111,err=111) res1 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock) c VAR: aa0thet is variable describing the average value of Foureir c VAR: expansion series c VAR: aathet is foureir expansion in theta/2 angle for full formula c VAR: look at the fitting equation in Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished read (ithep,*,end=111,err=111) &(aathet(l,i,j,k,iblock),l=1,ntheterm) read (ithep,*,end=111,err=111) & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle), & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle), & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle), & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle), & ll=1,ntheterm2) read (ithep,*,end=111,err=111) & (((ffthet(llll,lll,ll,i,j,k,iblock), & ffthet(lll,llll,ll,i,j,k,iblock), & ggthet(llll,lll,ll,i,j,k,iblock), & ggthet(lll,llll,ll,i,j,k,iblock), & 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 IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT C RECOMENTDED AFTER VERSION 3.3) c do i=1,nthetyp c do j=1,nthetyp c do l=1,ntheterm c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock) c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock) c enddo c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock) c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock) c enddo c do l=1,ntheterm c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock) c enddo c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock) c enddo c enddo C AND COMMENT THE LOOPS BELOW do i=1,nthetyp do j=1,nthetyp do l=1,ntheterm aathet(l,i,j,nthetyp+1,iblock)=0.0d0 aathet(l,nthetyp+1,i,j,iblock)=0.0d0 enddo aa0thet(i,j,nthetyp+1,iblock)=0.0d0 aa0thet(nthetyp+1,i,j,iblock)=0.0d0 enddo do l=1,ntheterm aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0 enddo aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0 enddo enddo C TILL HERE C Substitution for D aminoacids from symmetry. do iblock=1,2 do i=-nthetyp,0 do j=-nthetyp,nthetyp do k=-nthetyp,nthetyp aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock) do l=1,ntheterm aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock) enddo do ll=1,ntheterm2 do lll=1,nsingle bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock) ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock) ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock) eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock) enddo enddo do ll=1,ntheterm3 do lll=2,ndouble do llll=1,lll-1 ffthet(llll,lll,ll,i,j,k,iblock)= & ffthet(llll,lll,ll,-i,-j,-k,iblock) ffthet(lll,llll,ll,i,j,k,iblock)= & ffthet(lll,llll,ll,-i,-j,-k,iblock) ggthet(llll,lll,ll,i,j,k,iblock)= & -ggthet(llll,lll,ll,-i,-j,-k,iblock) ggthet(lll,llll,ll,i,j,k,iblock)= & -ggthet(lll,llll,ll,-i,-j,-k,iblock) enddo !ll enddo !lll enddo !llll enddo !k enddo !j enddo !i enddo !iblock 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 iblock=1,2 do i=0,nthetyp do j=-nthetyp,nthetyp do k=-nthetyp,nthetyp write (iout,'(//4a)') & 'Type ',toronelet(i),toronelet(j),toronelet(k) write (iout,'(//a,10x,a)') " l","a[l]" write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock) write (iout,'(i2,1pe15.5)') & (l,aathet(l,i,j,k,iblock),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,iblock),ccthet(m,l,i,j,k,iblock), & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock) 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,iblock), & ffthet(m,n,l,i,j,k,iblock), & ggthet(n,m,l,i,j,k,iblock), & ggthet(m,n,l,i,j,k,iblock) enddo enddo enddo enddo enddo enddo enddo call flush(iout) endif ELSE C here will be the apropriate recalibrating for D-aminoacid read (ithep,*,end=121,err=121) nthetyp do i=-nthetyp+1,nthetyp-1 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i) do j=0,nbend_kcc_Tb(i) read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i) enddo enddo if (lprint) then write (iout,'(a)') & "Parameters of the valence-only potentials" do i=-nthetyp+1,nthetyp-1 write (iout,'(2a)') "Type ",toronelet(i) do k=0,nbend_kcc_Tb(i) write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i) enddo enddo endif ENDIF ! TOR_MODE c write (2,*) "Start reading THETA_PDB",ithep_pdb do i=1,ntyp c write (2,*) 'i=',i 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 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 c 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) 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 write (2,*) "Start reading ROTAM_PDB" 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) write (2,*) "End reading ROTAM_PDB" #endif close(irotam) 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 SPLIT_FOURIERTOR = nloctyp.lt.0 nloctyp = iabs(nloctyp) #ifdef NEWCORR read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp) read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1) itype2loc(ntyp1)=nloctyp iloctyp(nloctyp)=ntyp1 do i=1,ntyp1 itype2loc(-i)=-itype2loc(i) enddo #else iloctyp(0)=10 iloctyp(1)=9 iloctyp(2)=20 iloctyp(3)=ntyp1 #endif do i=1,nloctyp iloctyp(-i)=-iloctyp(i) enddo c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1) c write (iout,*) "nloctyp",nloctyp, c & " iloctyp",(iloctyp(i),i=0,nloctyp) #ifdef NEWCORR bnew1=0.0d0 bnew2=0.0d0 ccnew=0.0d0 ddnew=0.0d0 eenew=0.0d0 e0new=0.0d0 do i=0,nloctyp-1 c write (iout,*) "NEWCORR",i read (ifourier,*,end=115,err=115) do ii=1,3 do j=1,2 read (ifourier,*,end=115,err=115) bnew1(ii,j,i) enddo enddo c write (iout,*) "NEWCORR BNEW1" c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2) do ii=1,3 do j=1,2 read (ifourier,*,end=115,err=115) bnew2(ii,j,i) enddo enddo c write (iout,*) "NEWCORR BNEW2" c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2) do kk=1,3 read (ifourier,*,end=115,err=115) ccnew(kk,1,i) read (ifourier,*,end=115,err=115) ccnew(kk,2,i) enddo c write (iout,*) "NEWCORR CCNEW" c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2) do kk=1,3 read (ifourier,*,end=115,err=115) ddnew(kk,1,i) read (ifourier,*,end=115,err=115) ddnew(kk,2,i) enddo c write (iout,*) "NEWCORR DDNEW" c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2) do ii=1,2 do jj=1,2 do kk=1,2 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i) enddo enddo enddo c write (iout,*) "NEWCORR EENEW1" c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2) do ii=1,3 read (ifourier,*,end=115,err=115) e0new(ii,i) enddo c write (iout,*) (e0new(ii,i),ii=1,3) enddo c write (iout,*) "NEWCORR EENEW" do i=0,nloctyp-1 do ii=1,3 ccnew(ii,1,i)=ccnew(ii,1,i)/2 ccnew(ii,2,i)=ccnew(ii,2,i)/2 ddnew(ii,1,i)=ddnew(ii,1,i)/2 ddnew(ii,2,i)=ddnew(ii,2,i)/2 enddo enddo do i=1,nloctyp-1 do ii=1,3 bnew1(ii,1,-i)= bnew1(ii,1,i) bnew1(ii,2,-i)=-bnew1(ii,2,i) bnew2(ii,1,-i)= bnew2(ii,1,i) bnew2(ii,2,-i)=-bnew2(ii,2,i) enddo do ii=1,3 c ccnew(ii,1,i)=ccnew(ii,1,i)/2 c ccnew(ii,2,i)=ccnew(ii,2,i)/2 c ddnew(ii,1,i)=ddnew(ii,1,i)/2 c ddnew(ii,2,i)=ddnew(ii,2,i)/2 ccnew(ii,1,-i)=ccnew(ii,1,i) ccnew(ii,2,-i)=-ccnew(ii,2,i) ddnew(ii,1,-i)=ddnew(ii,1,i) ddnew(ii,2,-i)=-ddnew(ii,2,i) enddo e0new(1,-i)= e0new(1,i) e0new(2,-i)=-e0new(2,i) e0new(3,-i)=-e0new(3,i) do kk=1,2 eenew(kk,1,1,-i)= eenew(kk,1,1,i) eenew(kk,1,2,-i)=-eenew(kk,1,2,i) eenew(kk,2,1,-i)=-eenew(kk,2,1,i) eenew(kk,2,2,-i)= eenew(kk,2,2,i) enddo enddo if (lprint) then write (iout,'(a)') "Coefficients of the multibody terms" c do i=-nloctyp+1,nloctyp-1 do i=-nloctyp,nloctyp write (iout,*) "Type: ",onelet(iloctyp(i)) write (iout,*) "Coefficients of the expansion of B1" do j=1,2 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of B2" do j=1,2 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of C" write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3) write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of D" write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3) write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of E" write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3) do j=1,2 do k=1,2 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2) enddo enddo enddo endif IF (SPLIT_FOURIERTOR) THEN do i=0,nloctyp-1 c write (iout,*) "NEWCORR TOR",i read (ifourier,*,end=115,err=115) do ii=1,3 do j=1,2 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i) enddo enddo c write (iout,*) "NEWCORR BNEW1 TOR" c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2) do ii=1,3 do j=1,2 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i) enddo enddo c write (iout,*) "NEWCORR BNEW2 TOR" c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2) do kk=1,3 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i) read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i) enddo c write (iout,*) "NEWCORR CCNEW TOR" c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2) do kk=1,3 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i) read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i) enddo c write (iout,*) "NEWCORR DDNEW TOR" c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2) do ii=1,2 do jj=1,2 do kk=1,2 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i) enddo enddo enddo c write (iout,*) "NEWCORR EENEW1 TOR" c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2) do ii=1,3 read (ifourier,*,end=115,err=115) e0newtor(ii,i) enddo c write (iout,*) (e0newtor(ii,i),ii=1,3) enddo c write (iout,*) "NEWCORR EENEW TOR" do i=0,nloctyp-1 do ii=1,3 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2 enddo enddo do i=1,nloctyp-1 do ii=1,3 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i) bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i) bnew2tor(ii,1,-i)= bnew2tor(ii,1,i) bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i) enddo do ii=1,3 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i) ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i) ddnewtor(ii,1,-i)=ddnewtor(ii,1,i) ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i) enddo e0newtor(1,-i)= e0newtor(1,i) e0newtor(2,-i)=-e0newtor(2,i) e0newtor(3,-i)=-e0newtor(3,i) do kk=1,2 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i) eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i) eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i) eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i) enddo enddo if (lprint) then write (iout,'(a)') & "Single-body coefficients of the torsional potentials" do i=-nloctyp+1,nloctyp-1 write (iout,*) "Type: ",onelet(iloctyp(i)) write (iout,*) "Coefficients of the expansion of B1tor" do j=1,2 write (iout,'(3hB1(,i1,1h),3f10.5)') & j,(bnew1tor(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of B2tor" do j=1,2 write (iout,'(3hB2(,i1,1h),3f10.5)') & j,(bnew2tor(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of Ctor" write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3) write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of Dtor" write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3) write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of Etor" write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3) do j=1,2 do k=1,2 write (iout,'(1hE,2i1,2f10.5)') & j,k,(eenewtor(l,j,k,i),l=1,2) enddo enddo enddo endif ELSE do i=-nloctyp+1,nloctyp-1 do ii=1,3 do j=1,2 bnew1tor(ii,j,i)=bnew1(ii,j,i) enddo enddo do ii=1,3 do j=1,2 bnew2tor(ii,j,i)=bnew2(ii,j,i) enddo enddo do ii=1,3 ccnewtor(ii,1,i)=ccnew(ii,1,i) ccnewtor(ii,2,i)=ccnew(ii,2,i) ddnewtor(ii,1,i)=ddnew(ii,1,i) ddnewtor(ii,2,i)=ddnew(ii,2,i) enddo enddo ENDIF !SPLIT_FOURIER_TOR #else if (lprint) & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)" do i=0,nloctyp-1 read (ifourier,*,end=115,err=115) read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13) if (lprint) then write (iout,*) 'Type ',onelet(iloctyp(i)) write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13) endif if (i.gt.0) then b(2,-i)= b(2,i) b(3,-i)= b(3,i) b(4,-i)=-b(4,i) b(5,-i)=-b(5,i) endif c B1(1,i) = b(3) c B1(2,i) = b(5) c B1(1,-i) = b(3) c B1(2,-i) = -b(5) c b1(1,i)=0.0d0 c b1(2,i)=0.0d0 c B1tilde(1,i) = b(3) c B1tilde(2,i) =-b(5) c B1tilde(1,-i) =-b(3) c B1tilde(2,-i) =b(5) c b1tilde(1,i)=0.0d0 c b1tilde(2,i)=0.0d0 c B2(1,i) = b(2) c B2(2,i) = b(4) c B2(1,-i) =b(2) c B2(2,-i) =-b(4) cc B1tilde(1,i) = b(3,i) cc B1tilde(2,i) =-b(5,i) C B1tilde(1,-i) =-b(3,i) C B1tilde(2,-i) =b(5,i) cc b1tilde(1,i)=0.0d0 cc b1tilde(2,i)=0.0d0 cc B2(1,i) = b(2,i) cc B2(2,i) = b(4,i) C B2(1,-i) =b(2,i) C B2(2,-i) =-b(4,i) c b2(1,i)=0.0d0 c b2(2,i)=0.0d0 CCold(1,1,i)= b(7,i) CCold(2,2,i)=-b(7,i) CCold(2,1,i)= b(9,i) CCold(1,2,i)= b(9,i) CCold(1,1,-i)= b(7,i) CCold(2,2,-i)=-b(7,i) CCold(2,1,-i)=-b(9,i) CCold(1,2,-i)=-b(9,i) 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 c Ctilde(1,1,i)= CCold(1,1,i) c Ctilde(1,2,i)= CCold(1,2,i) c Ctilde(2,1,i)=-CCold(2,1,i) c Ctilde(2,2,i)=-CCold(2,2,i) 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 DDold(1,1,i)= b(6,i) DDold(2,2,i)=-b(6,i) DDold(2,1,i)= b(8,i) DDold(1,2,i)= b(8,i) DDold(1,1,-i)= b(6,i) DDold(2,2,-i)=-b(6,i) DDold(2,1,-i)=-b(8,i) DDold(1,2,-i)=-b(8,i) 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 c Dtilde(1,1,i)= DD(1,1,i) c Dtilde(1,2,i)= DD(1,2,i) c Dtilde(2,1,i)=-DD(2,1,i) c Dtilde(2,2,i)=-DD(2,2,i) 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 EEold(1,1,i)= b(10,i)+b(11,i) EEold(2,2,i)=-b(10,i)+b(11,i) EEold(2,1,i)= b(12,i)-b(13,i) EEold(1,2,i)= b(12,i)+b(13,i) EEold(1,1,-i)= b(10,i)+b(11,i) EEold(2,2,-i)=-b(10,i)+b(11,i) EEold(2,1,-i)=-b(12,i)+b(13,i) EEold(1,2,-i)=-b(12,i)-b(13,i) write(iout,*) "TU DOCHODZE" print *,"JESTEM" 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 write (iout,*) write (iout,*) &"Coefficients of the cumulants (independent of valence angles)" do i=-nloctyp+1,nloctyp-1 write (iout,*) 'Type ',onelet(iloctyp(i)) write (iout,*) 'B1' write(iout,'(2f10.5)') B(3,i),B(5,i) write (iout,*) 'B2' write(iout,'(2f10.5)') B(2,i),B(4,i) write (iout,*) 'CC' do j=1,2 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i) enddo write(iout,*) 'DD' do j=1,2 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i) enddo write(iout,*) 'EE' do j=1,2 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i) enddo enddo endif #endif #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 IF (TOR_MODE.eq.0) THEN read (itorp,*,end=113,err=113) ntortyp read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp) do i=1,ntyp1 itype2loc(i)=itortyp(i) enddo do i=1,ntyp1 itype2loc(-i)=-itype2loc(i) enddo itortyp(ntyp1)=ntortyp do iblock=1,2 do i=-ntyp,-1 itortyp(i)=-itortyp(-i) enddo 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/)') 'Parameters of the SCCOR potentials:' do iblock=1,2 do i=0,ntortyp-1 do j=-ntortyp+1,ntortyp-1 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 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 Martix 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 Martix 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=0,ntortyp-1 do j=-ntortyp+1,ntortyp-1 do k=-ntortyp+1,ntortyp-1 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), & v1s(1,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 ELSE IF (TOR_MODE.eq.1) THEN C read valence-torsional parameters read (itorp,*,end=121,err=121) ntortyp nkcctyp=ntortyp write (iout,*) "Valence-torsional parameters read in ntortyp", & ntortyp read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp) write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp) #ifndef NEWCORR do i=1,ntyp1 itype2loc(i)=itortyp(i) enddo #endif do i=-ntyp,-1 itortyp(i)=-itortyp(-i) enddo do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 C first we read the cos and sin gamma parameters read (itorp,'(13x,a)',end=121,err=121) string write (iout,*) i,j,string read (itorp,*,end=121,err=121) & nterm_kcc(j,i),nterm_kcc_Tb(j,i) C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i) do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) read (itorp,*,end=121,err=121) ii,jj,kk, & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) enddo enddo enddo enddo enddo ELSE #ifdef NEWCORR c AL 4/8/16: Calculate coefficients from one-body parameters ntortyp=nloctyp do i=-ntyp1,ntyp1 itortyp(i)=itype2loc(i) enddo write (iout,*) &"Val-tor parameters calculated from cumulant coefficients ntortyp" & ,ntortyp do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 nterm_kcc(j,i)=2 nterm_kcc_Tb(j,i)=3 do k=1,nterm_kcc_Tb(j,i) do l=1,nterm_kcc_Tb(j,i) v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j) & +bnew1tor(k,2,i)*bnew2tor(l,2,j) v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j) & +bnew1tor(k,2,i)*bnew2tor(l,1,j) enddo enddo do k=1,nterm_kcc_Tb(j,i) do l=1,nterm_kcc_Tb(j,i) #ifdef CORRCD v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) & -ccnewtor(k,2,i)*ddnewtor(l,2,j)) v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) & +ccnewtor(k,1,i)*ddnewtor(l,2,j)) #else v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) & -ccnewtor(k,2,i)*ddnewtor(l,2,j)) v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) & +ccnewtor(k,1,i)*ddnewtor(l,2,j)) #endif enddo enddo c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma) enddo enddo #else write (iout,*) "TOR_MODE>1 only with NEWCORR" stop #endif ENDIF ! TOR_MODE if (tor_mode.gt.0 .and. lprint) then c Print valence-torsional parameters write (iout,'(a)') & "Parameters of the valence-torsional potentials" do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j) write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc" do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) write (iout,'(3i5,2f15.4)') & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) enddo enddo 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=119,err=119) nsccortyp #ifdef SCCORPDB read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp) do i=-ntyp,-1 isccortyp(i)=-isccortyp(-i) enddo iscprol=isccortyp(20) 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,*,end=119,err=119) &nterm_sccor(i,j),nlor_sccor(i,j) v0ijsccor=0.0d0 v0ijsccor1=0.0d0 v0ijsccor2=0.0d0 v0ijsccor3=0.0d0 si=-1.0d0 nterm_sccor(-i,j)=nterm_sccor(i,j) nterm_sccor(-i,-j)=nterm_sccor(i,j) nterm_sccor(i,-j)=nterm_sccor(i,j) do k=1,nterm_sccor(i,j) read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j) & ,v2sccor(k,l,i,j) if (j.eq.iscprol) then if (i.eq.isccortyp(10)) then v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) else v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 & +v2sccor(k,l,i,j)*dsqrt(0.75d0) v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 & +v1sccor(k,l,i,j)*dsqrt(0.75d0) v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) endif else if (i.eq.isccortyp(10)) then v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) else if (j.eq.isccortyp(10)) then v1sccor(k,l,-i,j)=v1sccor(k,l,i,j) v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j) else v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j) v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) endif endif endif v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j) v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j) v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j) si=-si enddo do k=1,nlor_sccor(i,j) read (isccor,*,end=119,err=119) 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(l,i,j)=v0ijsccor v0sccor(l,-i,j)=v0ijsccor1 v0sccor(l,i,-j)=v0ijsccor2 v0sccor(l,-i,-j)=v0ijsccor3 enddo enddo enddo close (isccor) #else read (isccor,*,end=119,err=119) (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,*,end=119,err=119) & nterm_sccor(i,j),nlor_sccor(i,j) v0ijsccor=0.0d0 si=-1.0d0 do k=1,nterm_sccor(i,j) read (isccor,*,end=119,err=119) 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,*,end=119,err=119) 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(l,i,j)=v0ijsccor enddo enddo enddo close (isccor) #endif 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 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 c lprint=.true. if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j), & ael6(i,j),ael3(i,j) c lprint=.false. 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=117,err=117)((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=117,err=117)((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 do i=1,ntyp read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp) enddo read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp) read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp) read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp) read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp) C now we start reading lipid do i=1,ntyp read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp) if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp) C print *,"WARNING!!" C do j=1,ntyp C epslip(i,j)=epslip(i,j)+0.05d0 C enddo enddo if (lprint) write(iout,*) epslip(1,1),"OK?" 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=117,err=117)((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) epslip(i,j)=epslip(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_aq(i,j)=epsij*rrij*rrij bb_aq(i,j)=-sigeps*epsij*rrij aa_aq(j,i)=aa_aq(i,j) bb_aq(j,i)=bb_aq(i,j) epsijlip=epslip(i,j) sigeps=dsign(1.0D0,epsijlip) epsijlip=dabs(epsijlip) aa_lip(i,j)=epsijlip*rrij*rrij bb_lip(i,j)=-sigeps*epsijlip*rrij aa_lip(j,i)=aa_lip(i,j) bb_lip(j,i)=bb_lip(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_aq(i,j),bb_aq(i,j),augm(i,j), & sigma(i,j),r0(i,j),chi(i,j),chi(j,i) endif enddo enddo #ifdef TUBE write(iout,*) "tube param" read(itube,*) epspeptube,sigmapeptube sigmapeptube=sigmapeptube**6 sigeps=dsign(1.0D0,epspeptube) epspeptube=dabs(epspeptube) pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep do i=1,ntyp read(itube,*) epssctube,sigmasctube sigmasctube=sigmasctube**6 sigeps=dsign(1.0D0,epssctube) epssctube=dabs(epssctube) sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i) enddo #endif #ifdef OLDSCP C C Define the SC-p interaction constants (hard-coded; old style) C do i=1,ntyp 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 c lprint=.true. if (lprint) then write (iout,*) "Parameters of SC-p interactions:" do i=1,ntyp 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 c lprint=.false. #endif C C Define the constants of the disulfide bridge C C ebr=-12.00D0 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 c akcm=0.0d0 c akth=0.0d0 c akct=0.0d0 c v1ss=0.0d0 c v2ss=0.0d0 c v3ss=0.0d0 C if(me.eq.king) then C write (iout,'(/a)') "Disulfide bridge parameters:" C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, C & ' v3ss:',v3ss C endif if (shield_mode.gt.0) then C VSolvSphere the volume of solving sphere C rpp(1,1) is the energy r0 for peptide group contact and will be used for it C there will be no distinction between proline peptide group and normal peptide C group in case of shielding parameters c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div", & VSolvSphere_div C long axis of side chain do i=1,ntyp long_r_sidechain(i)=vbldsc0(1,i) short_r_sidechain(i)=sigma0(i) enddo buff_shield=1.0d0 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 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip" 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" go to 999 121 write (iout,*) "Error in Czybyshev 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 C set the variables used for shielding effect C if (shield_mode.gt.0) then C VSolvSphere the volume of solving sphere C print *,pi,"pi" C rpp(1,1) is the energy r0 for peptide group contact and will be used for it C there will be no distinction between proline peptide group and normal peptide C group in case of shielding parameters C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 C long axis of side chain C do i=1,ntyp C long_r_sidechain(i)=vbldsc0(1,i) C short_r_sidechain(i)=sigma0(i) C enddo C lets set the buffor value C buff_shield=1.0d0 C endif return end c--------------------------------------------------------------------------- subroutine weightread 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 none 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' include 'COMMON.CONTROL' include 'COMMON.SHIELD' character*1000 weightcard integer i,j double precision scalscp,wlong c c READ energy-term weights c call card_concat(weightcard) call reada(weightcard,'WLONG',wlong,1.0D0) call reada(weightcard,'WSC',wsc,wlong) call reada(weightcard,'WSCP',wscp,wlong) call reada(weightcard,'WELEC',welec,1.0D0) call reada(weightcard,'WVDWPP',wvdwpp,welec) call reada(weightcard,'WEL_LOC',wel_loc,1.0D0) call reada(weightcard,'WCORR4',wcorr4,0.0D0) call reada(weightcard,'WCORR5',wcorr5,0.0D0) call reada(weightcard,'WCORR6',wcorr6,0.0D0) call reada(weightcard,'WTURN3',wturn3,1.0D0) call reada(weightcard,'WTURN4',wturn4,1.0D0) call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) call reada(weightcard,'WDFAD',wdfa_dist,0.0d0) call reada(weightcard,'WDFAT',wdfa_tor,0.0d0) call reada(weightcard,'WDFAN',wdfa_nei,0.0d0) call reada(weightcard,'WDFAB',wdfa_beta,0.0d0) call reada(weightcard,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) call reada(weightcard,'WSHIELD',wshield,1.0d0) call reada(weightcard,'WLT',wliptran,0.0D0) call reada(weightcard,'WTUBE',wtube,1.0D0) call reada(weightcard,'WSAXS',wsaxs,1.0D0) if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) if (wcorr4.gt.0.0d0) wcorr=wcorr4 weights(1)=wsc weights(2)=wscp weights(3)=welec weights(4)=wcorr weights(5)=wcorr5 weights(6)=wcorr6 weights(7)=wel_loc weights(8)=wturn3 weights(9)=wturn4 weights(10)=wturn6 weights(11)=wang weights(12)=wscloc weights(13)=wtor weights(14)=wtor_d weights(15)=wstrain weights(16)=wvdwpp weights(17)=wbond weights(18)=scal14 weights(21)=wsccor weights(26)=wsaxs weights(28)=wdfa_dist weights(29)=wdfa_tor weights(30)=wdfa_nei weights(31)=wdfa_beta do i=1,ntyp aad(i,1)=scalscp*aad(i,1) aad(i,2)=scalscp*aad(i,2) bad(i,1)=scalscp*bad(i,1) bad(i,2)=scalscp*bad(i,2) enddo if(me.eq.king.or..not.out1file) & write (iout,100) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, & wturn4,wturn6 100 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ & 'WELEC= ',f10.6,' (p-p electr)'/ & 'WVDWPP= ',f10.6,' (p-p VDW)'/ & 'WBOND= ',f10.6,' (stretching)'/ & 'WANG= ',f10.6,' (bending)'/ & 'WSCLOC= ',f10.6,' (SC local)'/ & 'WTOR= ',f10.6,' (torsional)'/ & 'WTORD= ',f10.6,' (double torsional)'/ & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ & 'WCORR4= ',f10.6,' (multi-body 4th order)'/ & 'WCORR5= ',f10.6,' (multi-body 5th order)'/ & 'WCORR6= ',f10.6,' (multi-body 6th order)'/ & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ & 'WTURN6= ',f10.6,' (turns, 6th order)') if(me.eq.king.or..not.out1file)then if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', & 'between contact pairs of peptide groups' write (iout,'(2(a,f5.3/))') & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr, & 'Range of quenching the correlation terms:',2*delt_corr else if (wcorr.gt.0.0d0) then write (iout,'(/2a/)') 'Hydrogen-bonding correlation ', & 'between contact pairs of peptide groups' endif write (iout,'(a,f8.3)') & 'Scaling factor of 1,4 SC-p interactions:',scal14 write (iout,'(a,f8.3)') & 'General scaling factor of SC-p interactions:',scalscp endif r0_corr=cutoff_corr-delt_corr wumb=1.0d0 call rescale_weights(t_bath) if(me.eq.king.or..not.out1file) & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc, #ifdef FOURBODY & wcorr,wcorr5,wcorr6, #endif & wsccor,wturn3, #ifdef FOURBODY & wturn4, #endif & wturn6 22 format (/'Energy-term weights (scaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ & 'WELEC= ',f10.6,' (p-p electr)'/ & 'WVDWPP= ',f10.6,' (p-p VDW)'/ & 'WBOND= ',f10.6,' (stretching)'/ & 'WANG= ',f10.6,' (bending)'/ & 'WSCLOC= ',f10.6,' (SC local)'/ & 'WTOR= ',f10.6,' (torsional)'/ & 'WTORD= ',f10.6,' (double torsional)'/ & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ #ifdef FOURBODY & 'WCORR4= ',f10.6,' (multi-body 4th order)'/ & 'WCORR5= ',f10.6,' (multi-body 5th order)'/ & 'WCORR6= ',f10.6,' (multi-body 6th order)'/ #endif & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ #ifdef FOURBODY & 'WTURN6= ',f10.6,' (turns, 6th order)' #endif & ) if(me.eq.king.or..not.out1file) & write (iout,*) "Reference temperature for weights calculation:", & temp0 #ifdef DFA write (iout,'(/a/)') "DFA pseudopotential parameters:" write (iout,'(a,f10.6,a)') & "WDFAD= ",wdfa_dist," (distance)", & "WDFAT= ",wdfa_tor," (backbone angles)", & "WDFAN= ",wdfa_nei," (neighbors)", & "WDFAB= ",wdfa_beta," (beta structure)" #endif call reada(weightcard,"D0CM",d0cm,3.78d0) call reada(weightcard,"AKCM",akcm,15.1d0) call reada(weightcard,"AKTH",akth,11.0d0) call reada(weightcard,"AKCT",akct,12.0d0) call reada(weightcard,"V1SS",v1ss,-1.08d0) call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) call reada(weightcard,"ATRISS",atriss,0.301D0) call reada(weightcard,"BTRISS",btriss,0.021D0) call reada(weightcard,"CTRISS",ctriss,1.001D0) call reada(weightcard,"DTRISS",dtriss,1.001D0) dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. enddo do i=1,max_cyst-1 do j=i+1,max_cyst dyn_ssbond_ij(i,j)=1.0d300 enddo enddo call flush(iout) call reada(weightcard,"HT",Ht,0.0D0) if (dyn_ss) then ss_depth=ebr/wsc-0.25*eps(1,1) Ht=Ht/wsc-0.25*eps(1,1) akcm=akcm*wstrain/wsc akth=akth*wstrain/wsc akct=akct*wstrain/wsc v1ss=v1ss*wstrain/wsc v2ss=v2ss*wstrain/wsc v3ss=v3ss*wstrain/wsc else if (wstrain.ne.0.0) then ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain else ss_depth=0.0 endif endif write (iout,*) "wshield,", wshield if(me.eq.king.or..not.out1file) then write (iout,*) "Parameters of the SS-bond potential:" write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth, & " AKCT",akct write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth write (iout,*)" HT",Ht write (iout,*) "Parameters of the 'trisulfide' potential" write (iout,*) "ATRISS=", atriss write (iout,*) "BTRISS=", btriss write (iout,*) "CTRISS=", ctriss write (iout,*) "DTRISS=", dtriss c print *,'indpdb=',indpdb,' pdbref=',pdbref endif return end