X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fparmread.F;h=ce1bf31546e0009f49871174002dc58676404b1b;hb=46f3f544d5a60812e2e11e03b619903f2ed052eb;hp=87be54f57237ad4498092bc12c941f1a160f2c3c;hpb=23c53ef4bfe4b994173cbbb61ecb859e829b4d48;p=unres.git diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F index 87be54f..ce1bf31 100644 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@ -21,15 +21,17 @@ C include 'COMMON.SCCOR' include 'COMMON.SCROT' include 'COMMON.FREE' + include 'COMMON.SHIELD' + include 'COMMON.CONTROL' 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 dimension blower(3,3,maxlob) - character*800 controlcard + character*900 controlcard character*256 bondname_t,thetname_t,rotname_t,torname_t, & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t, - & sccorname_t + & sccorname_t,tubename_t integer ilen external ilen character*16 key @@ -69,7 +71,13 @@ c If reading not own parameters, skip assignment call reada(controlcard,"ATRISS",atriss,0.3D0) call reada(controlcard,"BTRISS",btriss,0.02D0) call reada(controlcard,"CTRISS",ctriss,1.0D0) + call reada(controlcard,"LIPSCALE",lipscale,1.3D0) dyn_ss=(index(controlcard,'DYN_SS').gt.0) + write(iout,*) "ATRISS",atriss + write(iout,*) "BTRISS",btriss + write(iout,*) "CTRISS",ctriss + write(iout,*) "DTRISS",dtriss + C do i=1,maxres C dyn_ss_mask(i)=.false. C enddo @@ -136,6 +144,10 @@ c wbond=ww(18) wsccor=ww(19) whpb=ww(15) + wstrain=ww(15) + wliptran=ww(22) + wshield=ww(25) + wtube=ww(26) endif call card_concat(controlcard,.false.) @@ -174,6 +186,12 @@ c Return if not own parameters call reads(controlcard,"SCPPAR",scpname_t,scpname) open (iscpp,file=scpname_t,status='old') rewind(iscpp) + call reads(controlcard,"TUBEPAR",tubename_t,tubename) + write(iout,*) tubename_t + write(iout,*) tubename + open (itube,file=tubename_t,status='old') + rewind(itube) + write (iout,*) "Parameter set:",iparm write (iout,*) "Energy-term weights:" do i=1,n_ene @@ -207,7 +225,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 + read (ibond,*) vbldp0,vbldpdum,akp do i=1,ntyp nbondterm(i)=1 read (ibond,*) vbldsc0(1,i),aksc(1,i) @@ -219,7 +237,7 @@ c endif enddo #else - read (ibond,*) ijunk,vbldp0,akp,rjunk + read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk do i=1,ntyp read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i), & j=1,nbondterm(i)) @@ -245,6 +263,11 @@ c enddo enddo 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 @@ -389,7 +412,7 @@ C enddo enddo enddo - write (iout,*) "KURWA1" +C write (iout,*) "KURWA1" do iblock=1,2 do i=0,nthetyp do j=-nthetyp,nthetyp @@ -413,7 +436,7 @@ C enddo enddo enddo - write(iout,*) "KURWA1.1" +C write(iout,*) "KURWA1.1" C C For dummy ends assign glycine-type coefficients of theta-only terms; the C coefficients of theta-and-gamma-dependent terms are zero. @@ -433,7 +456,7 @@ C aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0 enddo enddo - write(iout,*) "KURWA1.5" +C write(iout,*) "KURWA1.5" C Substitution for D aminoacids from symmetry. do iblock=1,2 do i=-nthetyp,0 @@ -512,7 +535,7 @@ C call flush(iout) endif #endif - write(iout,*) 'KURWA2' +C write(iout,*) 'KURWA2' #ifdef CRYST_SC C C Read the parameters of the probability distribution/energy expression @@ -620,7 +643,7 @@ C enddo #endif close(irotam) - write (iout,*) 'KURWAKURWA' +C write (iout,*) 'KURWAKURWA' #ifdef CRYST_TOR C C Read torsional parameters in old format @@ -1095,6 +1118,14 @@ C---------------------- GB or BP potential ----------------------------- read (isidep,*)(sigii(i),i=1,ntyp) read (isidep,*)(chip(i),i=1,ntyp) read (isidep,*)(alp(i),i=1,ntyp) + do i=1,ntyp + read (isidep,*)(epslip(i,j),j=i,ntyp) +C write(iout,*) "WARNING!!",i,ntyp + write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp) +C do j=1,ntyp +C epslip(i,j)=epslip(i,j)+0.05d0 +C enddo + enddo C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@ -1133,6 +1164,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 @@ -1150,6 +1182,7 @@ C Calculate the "working" parameters of SC interactions. do i=1,ntyp do j=i,ntyp epsij=eps(i,j) + epsijlip=epslip(i,j) if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then rrij=sigma(i,j) else @@ -1161,10 +1194,16 @@ 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) + 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 @@ -1197,11 +1236,31 @@ 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 enddo + write(iout,*) "tube param" + read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, + & ccavtubpep,dcavtubpep,tubetranenepep + 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,acavtub(i),bcavtub(i), + & ccavtub(i),dcavtub(i),tubetranene(i) + 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 + C C Define the SC-p interaction constants C @@ -1274,7 +1333,7 @@ C V3SS = 13.7d0 write (iout,*) dyn_ss,'dyndyn' if (dyn_ss) then ss_depth=ebr/wsc-0.25*eps(1,1) - write(iout,*) akcm,whpb,wsc,'KURWA' +C write(iout,*) akcm,whpb,wsc,'KURWA' Ht=Ht/wsc-0.25*eps(1,1) akcm=akcm*whpb/wsc @@ -1295,5 +1354,22 @@ C if (lprint) then write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, & ' v3ss:',v3ss C endif + if (shield_mode.gt.0) then + pi=3.141592d0 +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 + 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_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 end