X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fparmread.F;h=43ceeac9097cd246aea442d9aa190d7bb193de30;hb=ad2687da4ea0c0ddde6c400a60c4817a689b3dff;hp=29c5dec2ef12a7084a72ae040e6e1b966f56be22;hpb=2951c9e609af806daf239c419202b721e55da3f5;p=unres.git diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 29c5dec..43ceeac 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -31,7 +31,7 @@ C character*1 toronelet(-2:2) /"p","a","G","A","P"/ logical lprint,LaTeX dimension blower(3,3,maxlob) - dimension b(13) +C dimension b(13) character*3 lancuch,ucase C C For printing parameters after they are read set the following in the UNRES @@ -59,7 +59,7 @@ c Read the virtual-bond parameters, masses, and moments of inertia c and Stokes' radii of the peptide group and side chains c #ifdef CRYST_BOND - read (ibond,*) vbldp0,akp,mp,ip,pstok + 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) @@ -71,7 +71,7 @@ c endif enddo #else - read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok + 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) @@ -243,12 +243,21 @@ C 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) @@ -269,7 +278,8 @@ C 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 +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 @@ -285,6 +295,7 @@ 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 @@ -300,7 +311,7 @@ c 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 @@ -377,6 +388,53 @@ C enddo call flush(iout) endif + 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 + write (2,*) "End reading THETA_PDB" + close (ithep_pdb) #endif close(ithep) #ifdef CRYST_SC @@ -470,7 +528,6 @@ 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 -CC TU JEST ZLE musibyc ntyp read (irotam,*,end=112,err=112) if (i.eq.10) then read (irotam,*,end=112,err=112) @@ -480,6 +537,50 @@ CC TU JEST ZLE musibyc ntyp 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) @@ -734,7 +835,7 @@ cc maxinter is maximum interaction sites si=-si enddo do k=1,nlor_sccor(i,j) - read (isccor,*,end=113,err=113) kk,vlor1sccor(k,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) @@ -748,26 +849,26 @@ cc maxinter is maximum interaction sites enddo close (isccor) #else - read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp) + 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=113,err=113) + 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=113,err=113) kk,v1sccor(k,l,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=113,err=113) kk,vlor1sccor(k,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) @@ -806,97 +907,106 @@ C write (iout,*) "Coefficients of the cumulants" endif read (ifourier,*) nloctyp + do i=0,nloctyp-1 read (ifourier,*,end=115,err=115) - read (ifourier,*,end=115,err=115) (b(ii),ii=1,13) + read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13) +#ifdef NEWCORR + read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3) + read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3) + read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1) + read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1) + read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1) +#endif if (lprint) then write (iout,*) 'Type',i - write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13) + write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13) endif - B1(1,i) = b(3) - B1(2,i) = b(5) - B1(1,-i) = b(3) - B1(2,-i) = -b(5) +c B1(1,i) = 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 - B1tilde(1,i) = b(3) - B1tilde(2,i) =-b(5) - B1tilde(1,-i) =-b(3) - B1tilde(2,-i) =b(5) +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 - B2(1,i) = b(2) - B2(2,i) = b(4) - B2(1,-i) =b(2) - B2(2,-i) =-b(4) +c B2(1,i) = b(2) +c B2(2,i) = b(4) +c B2(1,-i) =b(2) +c B2(2,-i) =-b(4) c b2(1,i)=0.0d0 c b2(2,i)=0.0d0 - CC(1,1,i)= b(7) - CC(2,2,i)=-b(7) - CC(2,1,i)= b(9) - CC(1,2,i)= b(9) - CC(1,1,-i)= b(7) - CC(2,2,-i)=-b(7) - CC(2,1,-i)=-b(9) - CC(1,2,-i)=-b(9) + 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) + 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) c CC(1,1,i)=0.0d0 c CC(2,2,i)=0.0d0 c CC(2,1,i)=0.0d0 c CC(1,2,i)=0.0d0 - Ctilde(1,1,i)=b(7) - Ctilde(1,2,i)=b(9) - Ctilde(2,1,i)=-b(9) - Ctilde(2,2,i)=b(7) - Ctilde(1,1,-i)=b(7) - Ctilde(1,2,-i)=-b(9) - Ctilde(2,1,-i)=b(9) - Ctilde(2,2,-i)=b(7) + 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) + 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) c Ctilde(1,1,i)=0.0d0 c Ctilde(1,2,i)=0.0d0 c Ctilde(2,1,i)=0.0d0 c Ctilde(2,2,i)=0.0d0 - DD(1,1,i)= b(6) - DD(2,2,i)=-b(6) - DD(2,1,i)= b(8) - DD(1,2,i)= b(8) - DD(1,1,-i)= b(6) - DD(2,2,-i)=-b(6) - DD(2,1,-i)=-b(8) - DD(1,2,-i)=-b(8) + 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) + 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) c DD(1,1,i)=0.0d0 c DD(2,2,i)=0.0d0 c DD(2,1,i)=0.0d0 c DD(1,2,i)=0.0d0 - Dtilde(1,1,i)=b(6) - Dtilde(1,2,i)=b(8) - Dtilde(2,1,i)=-b(8) - Dtilde(2,2,i)=b(6) - Dtilde(1,1,-i)=b(6) - Dtilde(1,2,-i)=-b(8) - Dtilde(2,1,-i)=b(8) - Dtilde(2,2,-i)=b(6) + 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) + 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) c Dtilde(1,1,i)=0.0d0 c Dtilde(1,2,i)=0.0d0 c Dtilde(2,1,i)=0.0d0 c Dtilde(2,2,i)=0.0d0 - EE(1,1,i)= b(10)+b(11) - EE(2,2,i)=-b(10)+b(11) - EE(2,1,i)= b(12)-b(13) - EE(1,2,i)= b(12)+b(13) - EE(1,1,-i)= b(10)+b(11) - EE(2,2,-i)=-b(10)+b(11) - EE(2,1,-i)=-b(12)+b(13) - EE(1,2,-i)=-b(12)-b(13) - + 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" 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 +c lprint=.true. if (lprint) then do i=1,nloctyp write (iout,*) 'Type',i @@ -914,10 +1024,11 @@ c ee(2,1,i)=ee(1,2,i) enddo write(iout,*) 'EE' do j=1,2 - write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i) + write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i) enddo enddo endif +c lprint=.false. C C Read electrostatic-interaction parameters @@ -940,8 +1051,10 @@ C 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 @@ -1140,7 +1253,7 @@ C 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 @@ -1148,6 +1261,7 @@ C & 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