X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fparmread.F;h=03efdb25b6b800a3004beac6efb4de7a2c01fee9;hb=39aaaf5f5ec4de449f18993aacdec4b43cabb1b0;hp=61b17cc021506e79de001a8f6fd4cc6050c84399;hpb=a12fec1b1bb1405b29ce7db86524648c70f5d635;p=unres.git diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 61b17cc..03efdb2 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -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,29 +907,37 @@ 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) +#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) 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 @@ -882,21 +991,22 @@ 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)+b(11) + EEold(2,2,i)=-b(10)+b(11) + EEold(2,1,i)= b(12)-b(13) + EEold(1,2,i)= b(12)+b(13) + EEold(1,1,-i)= b(10)+b(11) + EEold(2,2,-i)=-b(10)+b(11) + EEold(2,1,-i)=-b(12)+b(13) + EEold(1,2,-i)=-b(12)-b(13) + 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 - lprint=.true. +c lprint=.true. if (lprint) then write (iout,*) "Parameters of SC-p interactions:" do i=1,ntyp @@ -1148,7 +1261,7 @@ C & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2) enddo endif - lprint=.false. +c lprint=.false. #endif C C Define the constants of the disulfide bridge