X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_config.f90;fp=source%2Funres%2Fio_config.f90;h=0dce2e11063e1dd91b449b27860c8455841cf329;hb=5d299c1a16ab51f8206b8ee3b17c7bcabe9321b7;hp=451731d22c13615b35e9e7ebd01d85ea8bced6e0;hpb=b75e68533dc131221f1af10f6b2b81d65c1acade;p=unres4.git diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 451731d..0dce2e1 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -743,7 +743,7 @@ 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 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 +782,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 +800,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 +842,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 +1430,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 @@ -1916,7 +1928,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 +1992,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,10 +2039,13 @@ ! 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 + 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 @@ -2036,6 +2055,7 @@ 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 +2084,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,7 +2128,7 @@ 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 @@ -2523,11 +2551,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 +2620,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 +2775,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) @@ -2861,7 +2892,7 @@ character(len=640) :: controlcard real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,& - + integer i nglob_csa=0 eglob_csa=1d99 @@ -2900,6 +2931,14 @@ 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 call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) @@ -2907,6 +2946,25 @@ ! 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 +3066,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 +3743,8 @@ ! 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') ! print *,"Processor",myrank," opened file ISIDEP" ! print *,"Processor",myrank," opened parameter files" #elif (defined G77) @@ -3691,6 +3770,8 @@ 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') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& readonly) @@ -3721,6 +3802,8 @@ 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') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read')