X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_config.f90;h=f6a3919775187b62b4058af0d117bcc8a10b728f;hb=8ca97b16fe25b7053f258263899ba030572cc58f;hp=25a5e2438c011882975d46f48e071f86d90ef50e;hpb=bbbdc8e18680625d3004f414aec255e9ca7b7353;p=unres4.git diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 25a5e24..f6a3919 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -1412,6 +1412,24 @@ enddo endif enddo +!---------reading nucleic acid parameters for rotamers------------------- + allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp) + do i=1,ntyp_molec(2) + read (irotam_nucl,*,end=112,err=112) + do j=1,9 + read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i) + enddo + enddo + close(irotam_nucl) + if (lprint) then + write (iout,*) + write (iout,*) "Base rotamer parameters" + do i=1,ntyp_molec(2) + write (iout,'(a)') restyp(i,2) + write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9) + enddo + endif + ! ! Read the parameters of the probability distribution/energy expression ! of the side chains. @@ -1681,6 +1699,48 @@ enddo endif #endif + allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) + read (itorp_nucl,*,end=113,err=113) ntortyp_nucl +! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2) +!el from energy module--------- + allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2) + allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2) + + allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor) + allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) + allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor) + allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2) + + allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) + allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) +!el--------------------------- + nterm_nucl(:,:)=0 + nlor_nucl(:,:)=0 +!el-------------------- + read (itorp_nucl,*,end=113,err=113) & + (itortyp_nucl(i),i=1,ntyp_molec(2)) +! print *,itortyp_nucl(:) +!c write (iout,*) 'ntortyp',ntortyp + do i=1,ntortyp_nucl + do j=1,ntortyp_nucl + read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j) +! print *,nterm_nucl(i,j),nlor_nucl(i,j) + v0ij=0.0d0 + si=-1.0d0 + do k=1,nterm_nucl(i,j) + read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j) + v0ij=v0ij+si*v1_nucl(k,i,j) + si=-si + enddo + do k=1,nlor_nucl(i,j) + read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),& + vlor2_nucl(k,i,j),vlor3_nucl(k,i,j) + v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2) + enddo + v0_nucl(i,j)=v0ij + enddo + enddo + ! Read of Side-chain backbone correlation parameters ! Modified 11 May 2012 by Adasko !CC @@ -2013,6 +2073,7 @@ allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp) allocate(augm(ntyp,ntyp)) !(ntyp,ntyp) 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)) @@ -2122,6 +2183,7 @@ end select continue close (isidep) + !----------------------------------------------------------------------- ! Calculate the "working" parameters of SC interactions. @@ -2230,6 +2292,85 @@ endif enddo enddo + + allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2))) + allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2) + +! augm(:,:)=0.0D0 +! chip(:)=0.0D0 +! alp(:)=0.0D0 +! sigma0(:)=0.0D0 +! sigii(:)=0.0D0 +! rr0(:)=0.0D0 + + read (isidep_nucl,*) ipot_nucl +! print *,"TU?!",ipot_nucl + if (ipot_nucl.eq.1) then + do i=1,ntyp_molec(2) + do j=i,ntyp_molec(2) + read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),& + elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j) + enddo + enddo + else + do i=1,ntyp_molec(2) + do j=i,ntyp_molec(2) + read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),& + chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),& + elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j) + enddo + enddo + endif +! rpp(1,1)=2**(1.0/6.0)*5.16158 + do i=1,ntyp_molec(2) + do j=i,ntyp_molec(2) + rrij=sigma_nucl(i,j) + r0_nucl(i,j)=rrij + r0_nucl(j,i)=rrij + rrij=rrij**expon + epsij=4*eps_nucl(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_nucl(i,j)=epsij*rrij*rrij + bb_nucl(i,j)=-sigeps*epsij*rrij + ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij) + ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij + ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij + ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij + sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- & + 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i))) + enddo + do j=1,i-1 + aa_nucl(i,j)=aa_nucl(j,i) + bb_nucl(i,j)=bb_nucl(j,i) + ael3_nucl(i,j)=ael3_nucl(j,i) + ael6_nucl(i,j)=ael6_nucl(j,i) + ael63_nucl(i,j)=ael63_nucl(j,i) + ael32_nucl(i,j)=ael32_nucl(j,i) + elpp3_nucl(i,j)=elpp3_nucl(j,i) + elpp6_nucl(i,j)=elpp6_nucl(j,i) + elpp63_nucl(i,j)=elpp63_nucl(j,i) + elpp32_nucl(i,j)=elpp32_nucl(j,i) + eps_nucl(i,j)=eps_nucl(j,i) + sigma_nucl(i,j)=sigma_nucl(j,i) + sigmaii_nucl(i,j)=sigmaii_nucl(j,i) + enddo + enddo + write(iout,*) "tube param" read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, & ccavtubpep,dcavtubpep,tubetranenepep @@ -2302,6 +2443,20 @@ endif ! lprint=.false. #endif + allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2) + + do i=1,ntyp_molec(2) + read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i) + enddo + do i=1,ntyp_molec(2) + aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12 + bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6 + enddo + r0pp=1.12246204830937298142*5.16158 + epspp=4.95713/4 + AEES=108.661 + BEES=0.433246 + ! ! Define the constants of the disulfide bridge ! @@ -2396,7 +2551,7 @@ ! Read the PDB file and convert the peptide geometry into virtual-chain ! geometry. use geometry_data - use energy_data, only: itype + use energy_data, only: itype,istype use control_data use compare_data use MPI_data @@ -2433,6 +2588,7 @@ real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm) real(kind=8),dimension(:,:), allocatable :: c_temporary integer,dimension(:,:) , allocatable :: itype_temporary + integer,dimension(:),allocatable :: istype_temp efree_temp=0.0d0 ibeg=1 ishift1=0 @@ -2544,7 +2700,7 @@ ishift1=ishift1-1 !!!!! ! write (iout,*) "New chain started",ires,ishift,ishift1,"!" ires=ires-ishift+ishift1 - print *,ires,ishift,ishift1 +! print *,ires,ishift,ishift1 ires_old=ires ibeg=0 else @@ -2564,6 +2720,15 @@ molecule=2 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule) ! nres_molec(molecule)=nres_molec(molecule)+1 + read (card(19:19),'(a1)') sugar + isugar=sugarcode(sugar,ires) +! if (ibeg.eq.1) then +! istype(1)=isugar +! else + istype(ires)=isugar +! print *,"ires=",ires,istype(ires) +! endif + endif endif else @@ -2593,7 +2758,7 @@ .or. atom.eq."C4'" .or. atom.eq."O4'")) then read(card(31:54),'(3f8.3)') (ccc(j),j=1,3) !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3) - print *,ires,ishift,ishift1 +! print *,ires,ishift,ishift1 counter=counter+1 ! iii=iii+1 ! do j=1,3 @@ -2602,6 +2767,7 @@ do j=1,3 c(j,ires)=c(j,ires)+ccc(j)/5.0 enddo + print *,counter,molecule if (counter.eq.5) then ! iii=iii+1 nres_molec(molecule)=nres_molec(molecule)+1 @@ -2611,13 +2777,14 @@ counter=0 endif ! print *, "ATOM",atom(1:3) - else if (atom(1:3).eq."C5'") then + else if (atom.eq."C5'") then read (card(19:19),'(a1)') sugar isugar=sugarcode(sugar,ires) if (ibeg.eq.1) then istype(1)=isugar else istype(ires)=isugar +! print *,ires,istype(ires) endif if (unres_pdb) then read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) @@ -2683,7 +2850,7 @@ if (itype(i,k).eq.ntyp1_molec(k)) then if (itype(i+1,k).eq.ntyp1_molec(k)) then if (itype(i+2,k).eq.0) then - print *,"masz sieczke" +! print *,"masz sieczke" do j=1,5 if (itype(i+2,j).ne.0) then itype(i+1,k)=0 @@ -2707,7 +2874,7 @@ e2(2)=1.0d0 e2(3)=0.0d0 endif !fail - print *,i,'a tu?' +! print *,i,'a tu?' do j=1,3 c(j,i)=c(j,i-1)-1.9d0*e2(j) enddo @@ -2905,6 +3072,7 @@ allocate(c_temporary(3,2*nres)) allocate(itype_temporary(nres,5)) allocate(molnum(nres)) + allocate(istype_temp(nres)) itype_temporary(:,:)=0 seqalingbegin=1 do k=1,5 @@ -2916,6 +3084,7 @@ enddo itype_temporary(seqalingbegin,k)=itype(i,k) + istype_temp(seqalingbegin)=istype(i) molnum(seqalingbegin)=k seqalingbegin=seqalingbegin+1 endif @@ -2929,6 +3098,7 @@ do k=1,5 do i=1,nres itype(i,k)=itype_temporary(i,k) + istype(i)=istype_temp(i) enddo enddo if (itype(1,1).eq.ntyp1) then @@ -2966,7 +3136,7 @@ enddo endif - print *,seqalingbegin,nres +! print *,seqalingbegin,nres if(.not.allocated(vbld)) then allocate(vbld(2*nres)) do i=1,2*nres @@ -3045,6 +3215,7 @@ kkk=1 lll=0 cou=1 + write (iout,*) "symetr", symetr do i=1,nres lll=lll+1 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) @@ -3537,7 +3708,7 @@ ! character(len=80) :: ucase character(len=320) :: controlcard !el local variables - integer :: i + integer :: i,j real(kind=8) :: eta call card_concat(controlcard,.true.) @@ -3650,21 +3821,27 @@ if(me.eq.king.or..not.out1file) & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",& eta - gamp=scal_fric*(pstok(i)+rwat)*eta - stdfp=dsqrt(2*Rb*t_bath/d_time) - allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1) +! allocate(gamp + do j=1,5 !types of molecules + gamp(j)=scal_fric*(pstok(j)+rwat)*eta + stdfp(j)=dsqrt(2*Rb*t_bath/d_time) + enddo + allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1) + do j=1,5 !types of molecules do i=1,ntyp - gamsc(i)=scal_fric*(restok(i,1)+rwat)*eta - stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) + gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta + stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time) enddo + enddo + if(me.eq.king.or..not.out1file)then write (iout,'(/2a/)') & "Radii of site types and friction coefficients and std's of",& " stochastic forces of fully exposed sites" - write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp) + write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1)) do i=1,ntyp write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),& - gamsc(i),stdfsc(i)*dsqrt(gamsc(i)) + gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1)) enddo endif else if (tbf) then @@ -4205,6 +4382,17 @@ open (irotam_pdb,file=rotname_pdb,status='old',action='read') #endif #endif + call getenv_loc('SCPPAR_NUCL',scpname_nucl) +#if defined(WINIFL) || defined(WINPGI) + open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared) +#elif (defined CRAY) || (defined AIX) + open (iscpp_nucl,file=scpname_nucl,status='old',action='read') +#elif (defined G77) + open (iscpp_nucl,file=scpname_nucl,status='old') +#else + open (iscpp_nucl,file=scpname_nucl,status='old',action='read') +#endif + #ifndef OLDSCP ! ! 8/9/01 In the newest version SCp interaction constants are read from a file