X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fparmread.F;h=ae4d7100649e77a67a91f6cbbbd1a447c41b7bbc;hb=25618f9f83673a7063414fe1e17415d138f58da8;hp=03efdb25b6b800a3004beac6efb4de7a2c01fee9;hpb=0195ef1798fae1a859d7ef35fc619308bfae5dba;p=unres.git diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 03efdb2..ae4d710 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -26,12 +26,14 @@ C include 'COMMON.SBRIDGE' include 'COMMON.MD' 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 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 +61,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,8 +73,9 @@ 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 + print *,i 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) @@ -97,6 +100,12 @@ c enddo enddo endif +C reading lipid parameters + 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 @@ -910,7 +919,7 @@ C 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) @@ -920,7 +929,7 @@ C #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 c B1(1,i) = b(3) c B1(2,i) = b(5) @@ -938,68 +947,79 @@ c B2(1,i) = b(2) c B2(2,i) = b(4) c B2(1,-i) =b(2) c B2(2,-i) =-b(4) + B1tilde(1,i) = b(3,i) + B1tilde(2,i) =-b(5,i) +C B1tilde(1,-i) =-b(3,i) +C B1tilde(2,-i) =b(5,i) + b1tilde(1,i)=0.0d0 + b1tilde(2,i)=0.0d0 + B2(1,i) = b(2,i) + 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 - 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 - 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) + 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 @@ -1075,7 +1095,7 @@ C & ', exponents are ',expon,2*expon goto (10,20,30,30,40) ipot C----------------------- LJ potential --------------------------------- - 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), + 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:' @@ -1087,7 +1107,7 @@ C----------------------- LJ potential --------------------------------- endif goto 50 C----------------------- LJK potential -------------------------------- - 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), + 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:' @@ -1101,12 +1121,21 @@ C----------------------- LJK potential -------------------------------- goto 50 C---------------------- GB or BP potential ----------------------------- 30 do i=1,ntyp - read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp) + read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp) enddo read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp) read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp) read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp) read (isidep,*,end=116,err=116)(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) +C print *,"WARNING!!" +C do j=1,ntyp +C epslip(i,j)=epslip(i,j)+0.05d0 +C enddo + enddo + 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 @@ -1125,7 +1154,7 @@ C For the GB potential convert sigma'**2 into chi' endif goto 50 C--------------------- GBV potential ----------------------------------- - 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp), + 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 @@ -1145,6 +1174,7 @@ 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 @@ -1173,10 +1203,17 @@ C Calculate the "working" parameters of SC interactions. 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) + 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 @@ -1209,7 +1246,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) 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), + & 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 @@ -1266,7 +1303,7 @@ c lprint=.false. C C Define the constants of the disulfide bridge C - ebr=-5.50D0 +C ebr=-12.00D0 c c Old arbitrary potential - commented out. c @@ -1277,13 +1314,13 @@ c Constants of the disulfide-bond potential determined based on the RHF/6-31G** c energy surface of diethyl disulfide. c A. Liwo and U. Kozlowska, 11/24/03 c - D0CM = 3.78d0 - AKCM = 15.1d0 - AKTH = 11.0d0 - AKCT = 12.0d0 - V1SS =-1.08d0 - V2SS = 7.61d0 - V3SS = 13.7d0 +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 @@ -1291,14 +1328,33 @@ c v1ss=0.0d0 c v2ss=0.0d0 c v3ss=0.0d0 - if(me.eq.king) then - write (iout,'(/a)') "Disulfide bridge parameters:" - write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr - write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm - write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct - write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, - & ' v3ss:',v3ss - endif +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 +C set the variables used for shielding effect +C write (iout,*) "SHIELD MODE",shield_mode +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 write (iout,*) VSolvSphere,VSolvSphere_div +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 111 write (iout,*) "Error reading bending energy parameters." goto 999 @@ -1313,6 +1369,8 @@ c v3ss=0.0d0 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." @@ -1368,6 +1426,22 @@ 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