X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_config.f90;h=0c2d645a05b9409c40fc9dd34fc420e58a816889;hb=refs%2Fheads%2Fdevel;hp=451731d22c13615b35e9e7ebd01d85ea8bced6e0;hpb=b75e68533dc131221f1af10f6b2b81d65c1acade;p=unres4.git diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 451731d..0c2d645 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -30,7 +30,7 @@ ! !----------------------------------------------------------------------------- contains -#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) +#if !defined(WHAM_RUN) && !defined(CLUSTER) !----------------------------------------------------------------------------- ! bank.F io_csa !----------------------------------------------------------------------------- @@ -735,7 +735,7 @@ character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/) logical :: lprint,LaTeX real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob) - real(kind=8),dimension(13) :: b + real(kind=8),dimension(13) :: bN character(len=3) :: lancuch !,ucase !el local variables integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm @@ -743,7 +743,8 @@ real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,& dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,& sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,& - res1 + res1,epsijlip,epspeptube,epssctube,sigmapeptube, & + sigmasctube integer :: ichir1,ichir2 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)) @@ -782,7 +783,8 @@ allocate(isc(ntyp+1)) !(ntyp+1) allocate(restok(ntyp+1)) !(ntyp+1) allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp) - + allocate(long_r_sidechain(ntyp)) + allocate(short_r_sidechain(ntyp)) dsc(:)=0.0d0 dsc_inv(:)=0.0d0 @@ -799,7 +801,7 @@ endif enddo #else - read (ibond,*) junk,vbldpDUM,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) @@ -841,6 +843,17 @@ theta0(:)=0.0D0 sig0(:)=0.0D0 sigc0(:)=0.0D0 + allocate(liptranene(ntyp)) +!C reading lipid parameters + write (iout,*) "iliptranpar",iliptranpar + call flush(iout) + read(iliptranpar,*) pepliptran + print *,pepliptran + do i=1,ntyp + read(iliptranpar,*) liptranene(i) + print *,liptranene(i) + enddo + close(iliptranpar) #ifdef CRYST_THETA ! @@ -1418,8 +1431,8 @@ do i=-ntyp,-1 itortyp(i)=-itortyp(-i) enddo -! itortyp(ntyp1)=ntortyp -! itortyp(-ntyp1)=-ntortyp + itortyp(ntyp1)=ntortyp + itortyp(-ntyp1)=-ntortyp do iblock=1,2 write (iout,*) 'ntortyp',ntortyp do i=0,ntortyp-1 @@ -1769,88 +1782,88 @@ 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) (bN(ii),ii=1,13) 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)') ('bN(',ii,')=',bN(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) + B1(1,i) = bN(3) + B1(2,i) = bN(5) + B1(1,-i) = bN(3) + B1(2,-i) = -bN(5) ! b1(1,i)=0.0d0 ! 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) + B1tilde(1,i) = bN(3) + B1tilde(2,i) =-bN(5) + B1tilde(1,-i) =-bN(3) + B1tilde(2,-i) =bN(5) ! b1tilde(1,i)=0.0d0 ! 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) + B2(1,i) = bN(2) + B2(2,i) = bN(4) + B2(1,-i) =bN(2) + B2(2,-i) =-bN(4) ! b2(1,i)=0.0d0 ! 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)= bN(7) + CC(2,2,i)=-bN(7) + CC(2,1,i)= bN(9) + CC(1,2,i)= bN(9) + CC(1,1,-i)= bN(7) + CC(2,2,-i)=-bN(7) + CC(2,1,-i)=-bN(9) + CC(1,2,-i)=-bN(9) ! CC(1,1,i)=0.0d0 ! CC(2,2,i)=0.0d0 ! CC(2,1,i)=0.0d0 ! 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)=bN(7) + Ctilde(1,2,i)=bN(9) + Ctilde(2,1,i)=-bN(9) + Ctilde(2,2,i)=bN(7) + Ctilde(1,1,-i)=bN(7) + Ctilde(1,2,-i)=-bN(9) + Ctilde(2,1,-i)=bN(9) + Ctilde(2,2,-i)=bN(7) ! Ctilde(1,1,i)=0.0d0 ! Ctilde(1,2,i)=0.0d0 ! Ctilde(2,1,i)=0.0d0 ! 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)= bN(6) + DD(2,2,i)=-bN(6) + DD(2,1,i)= bN(8) + DD(1,2,i)= bN(8) + DD(1,1,-i)= bN(6) + DD(2,2,-i)=-bN(6) + DD(2,1,-i)=-bN(8) + DD(1,2,-i)=-bN(8) ! DD(1,1,i)=0.0d0 ! DD(2,2,i)=0.0d0 ! DD(2,1,i)=0.0d0 ! 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)=bN(6) + Dtilde(1,2,i)=bN(8) + Dtilde(2,1,i)=-bN(8) + Dtilde(2,2,i)=bN(6) + Dtilde(1,1,-i)=bN(6) + Dtilde(1,2,-i)=-bN(8) + Dtilde(2,1,-i)=bN(8) + Dtilde(2,2,-i)=bN(6) ! Dtilde(1,1,i)=0.0d0 ! Dtilde(1,2,i)=0.0d0 ! Dtilde(2,1,i)=0.0d0 ! 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) + EE(1,1,i)= bN(10)+bN(11) + EE(2,2,i)=-bN(10)+bN(11) + EE(2,1,i)= bN(12)-bN(13) + EE(1,2,i)= bN(12)+bN(13) + EE(1,1,-i)= bN(10)+bN(11) + EE(2,2,-i)=-bN(10)+bN(11) + EE(2,1,-i)=-bN(12)+bN(13) + EE(1,2,-i)=-bN(12)-bN(13) ! ee(1,1,i)=1.0d0 ! ee(2,2,i)=1.0d0 @@ -1916,7 +1929,7 @@ allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2) allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp) allocate(chip(ntyp1),alp(ntyp1)) !(ntyp) - + allocate(epslip(ntyp,ntyp)) augm(:,:)=0.0D0 chip(:)=0.0D0 alp(:)=0.0D0 @@ -1980,6 +1993,10 @@ 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) + do i=1,ntyp + read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp) + enddo + ! For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@ -2023,19 +2040,33 @@ ! Calculate the "working" parameters of SC interactions. !el from module energy - COMMON.INTERACT------- - allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) + allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) + allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp) allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1) - aa(:,:)=0.0D0 - bb(:,:)=0.0D0 + allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),& + dcavtub(ntyp1)) + allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),& + tubetranene(ntyp1)) + aa_aq(:,:)=0.0D0 + bb_aq(:,:)=0.0D0 + aa_lip(:,:)=0.0D0 + bb_lip(:,:)=0.0D0 chi(:,:)=0.0D0 sigma(:,:)=0.0D0 r0(:,:)=0.0D0 - + acavtub(:)=0.0d0 + bcavtub(:)=0.0d0 + ccavtub(:)=0.0d0 + dcavtub(:)=0.0d0 + sc_aa_tube_par(:)=0.0d0 + sc_bb_tube_par(:)=0.0d0 + !-------------------------------- 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 @@ -2064,10 +2095,18 @@ 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) +!C write(iout,*) aa_lip if (ipot.gt.2) then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 @@ -2100,11 +2139,30 @@ 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 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2) bad(:,:)=0.0D0 @@ -2523,11 +2581,11 @@ e2(3)=0.0d0 endif do j=1,3 - c(j,nres)=c(j,nres-1)-3.8d0*e2(j) + c(j,nres)=c(j,nres-1)-1.9d0*e2(j) enddo else do j=1,3 - dcj=c(j,nres-2)-c(j,nres-3) + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 c(j,nres)=c(j,nres-1)+dcj c(j,2*nres)=c(j,nres) enddo @@ -2592,11 +2650,11 @@ e2(3)=0.0d0 endif do j=1,3 - c(j,1)=c(j,2)-3.8d0*e2(j) + c(j,1)=c(j,2)-1.9d0*e2(j) enddo else do j=1,3 - dcj=c(j,4)-c(j,3) + dcj=(c(j,4)-c(j,3))/2.0 c(j,1)=c(j,2)-dcj c(j,nres+1)=c(j,1) enddo @@ -2747,6 +2805,9 @@ ! enddiagnostic ! makes copy of chains write (iout,*) "symetr", symetr + do j=1,3 + dc(j,0)=c(j,1) + enddo if (symetr.gt.1) then call permut(symetr) @@ -2811,7 +2872,7 @@ return end subroutine readpdb -#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) +#if !defined(WHAM_RUN) && !defined(CLUSTER) !----------------------------------------------------------------------------- ! readrtns_CSA.F !----------------------------------------------------------------------------- @@ -2861,7 +2922,7 @@ character(len=640) :: controlcard real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,& - + integer i nglob_csa=0 eglob_csa=1d99 @@ -2900,13 +2961,60 @@ timem=timlim modecalc=0 call reada(controlcard,"T_BATH",t_bath,300.0d0) +!C SHIELD keyword sets if the shielding effect of side-chains is used +!C 0 denots no shielding is used all peptide are equally despite the +!C solvent accesible area +!C 1 the newly introduced function +!C 2 reseved for further possible developement + call readi(controlcard,'SHIELD',shield_mode,0) +!C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write(iout,*) "shield_mode",shield_mode !C Varibles set size of box + with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 + write (iout,*) "with_theta_constr ",with_theta_constr + AFMlog=(index(controlcard,'AFM')) + selfguide=(index(controlcard,'SELFGUIDE')) + print *,'AFMlog',AFMlog,selfguide,"KUPA" + call readi(controlcard,'GENCONSTR',genconstr,0) call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) + call readi(controlcard,'TUBEMOD',tubemode,0) + write (iout,*) TUBEmode,"TUBEMODE" + if (TUBEmode.gt.0) then + call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) + call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) + call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0) + call reada(controlcard,"RTUBE",tubeR0,0.0d0) + call reada(controlcard,"TUBETOP",bordtubetop,boxzsize) + call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0) + call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0) + buftubebot=bordtubebot+tubebufthick + buftubetop=bordtubetop-tubebufthick + endif + ! CUTOFFF ON ELECTROSTATICS call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0) call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0) + write(iout,*) "R_CUT_ELE=",r_cut_ele +! Lipidic parameters + call reada(controlcard,"LIPTHICK",lipthick,0.0d0) + call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) + if (lipthick.gt.0.0d0) then + bordliptop=(boxzsize+lipthick)/2.0 + bordlipbot=bordliptop-lipthick + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & + write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" + buflipbot=bordlipbot+lipbufthick + bufliptop=bordliptop-lipbufthick + if ((lipbufthick*2.0d0).gt.lipthick) & + write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" + endif !lipthick.gt.0 + write(iout,*) "bordliptop=",bordliptop + write(iout,*) "bordlipbot=",bordlipbot + write(iout,*) "bufliptop=",bufliptop + write(iout,*) "buflipbot=",buflipbot + write (iout,*) "SHIELD MODE",shield_mode !C------------------------- minim=(index(controlcard,'MINIMIZE').gt.0) @@ -3008,6 +3116,25 @@ if(me.eq.king.or..not.out1file) & write (iout,'(2a)') diagmeth(kdiag),& ' routine used to diagonalize matrices.' + 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) + write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& + sigma0(i) + enddo + buff_shield=1.0d0 + endif return end subroutine read_control !----------------------------------------------------------------------------- @@ -3666,6 +3793,11 @@ ! print *,"Processor",myrank," opened file IELEP" call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',action='read') + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',action='read') + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old',action='read') + ! print *,"Processor",myrank," opened file ISIDEP" ! print *,"Processor",myrank," opened parameter files" #elif (defined G77) @@ -3691,6 +3823,10 @@ open (ielep,file=elename,status='old') call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old') + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old') + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& readonly) @@ -3721,6 +3857,11 @@ open (ielep,file=elename,status='old',readonly) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly) + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',action='read') + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old',action='read') + #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') @@ -3917,6 +4058,7 @@ open(irest2,file=rest2name,status='unknown') read(irest2,*) totT,EK,potE,totE,t_bath + totTafm=totT ! do i=1,2*nres ! AL 4/17/17: Now reading d_t(0,:) too do i=0,2*nres