X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_config.f90;h=f6a3919775187b62b4058af0d117bcc8a10b728f;hb=8ca97b16fe25b7053f258263899ba030572cc58f;hp=277b6ba44e4d091a71d90262c2e3993505d49515;hpb=834a82f8158d4d91c1890efe979bf66cc399aea4;p=unres4.git diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 277b6ba..f6a3919 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -776,12 +776,15 @@ ! allocate(dsc(ntyp1)) !(ntyp1) allocate(dsc_inv(ntyp1)) !(ntyp1) + allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp) + allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp) + allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp) allocate(nbondterm(ntyp)) !(ntyp) allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp) - allocate(msc(ntyp+1)) !(ntyp+1) - allocate(isc(ntyp+1)) !(ntyp+1) - allocate(restok(ntyp+1)) !(ntyp+1) + allocate(msc(ntyp+1,5)) !(ntyp+1) + allocate(isc(ntyp+1,5)) !(ntyp+1) + allocate(restok(ntyp+1,5)) !(ntyp+1) allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(long_r_sidechain(ntyp)) allocate(short_r_sidechain(ntyp)) @@ -801,10 +804,10 @@ endif enddo #else - read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok - do i=1,ntyp + read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1) + do i=1,ntyp_molec(1) read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),& - j=1,nbondterm(i)),msc(i),isc(i),restok(i) + j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1) dsc(i) = vbldsc0(1,i) if (i.eq.10) then dsc_inv(i)=0.0D0 @@ -813,14 +816,38 @@ endif enddo #endif + read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2) + do i=1,ntyp_molec(2) + nbondterm_nucl(i)=1 + read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2) +! dsc(i) = vbldsc0_nucl(1,i) +! if (i.eq.10) then +! dsc_inv(i)=0.0D0 +! else +! dsc_inv(i)=1.0D0/dsc(i) +! endif + enddo +! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2) +! do i=1,ntyp_molec(2) +! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),& +! aksc_nucl(j,i),abond0_nucl(j,i),& +! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2) +! dsc(i) = vbldsc0(1,i) +! if (i.eq.10) then +! dsc_inv(i)=0.0D0 +! else +! dsc_inv(i)=1.0D0/dsc(i) +! endif +! enddo + if (lprint) then write(iout,'(/a/)')"Dynamic constants of the interaction sites:" write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',& 'inertia','Pstok' - write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok + write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1) do i=1,ntyp write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),& - vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i) + vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1) do j=2,nbondterm(i) write (iout,'(13x,3f10.5)') & vbldsc0(j,i),aksc(j,i),abond0(j,i) @@ -1207,6 +1234,65 @@ close (ithep_pdb) #endif close(ithep) +!--------------- Reading theta parameters for nucleic acid------- + read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,& + ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl + nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl) + allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) + allocate(aa0thet_nucl(nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) +!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) +!(maxtheterm,-maxthetyp1:maxthetyp1,& +! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) +!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,& +! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1)) + +!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,& +! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) + + read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2)) + + aa0thet_nucl(:,:,:)=0.0d0 + aathet_nucl(:,:,:,:)=0.0d0 + bbthet_nucl(:,:,:,:,:)=0.0d0 + ccthet_nucl(:,:,:,:,:)=0.0d0 + ddthet_nucl(:,:,:,:,:)=0.0d0 + eethet_nucl(:,:,:,:,:)=0.0d0 + ffthet_nucl(:,:,:,:,:,:)=0.0d0 + ggthet_nucl(:,:,:,:,:,:)=0.0d0 + + do i=1,nthetyp_nucl + do j=1,nthetyp_nucl + do k=1,nthetyp_nucl + read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3 + read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k) + read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl) + read (ithep_nucl,*,end=111,err=111) & + (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & + (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & + (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & + (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl) + read (ithep_nucl,*,end=111,err=111) & + (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), & + ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), & + llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl) + enddo + enddo + enddo !------------------------------------------- allocate(nlob(ntyp1)) !(ntyp1) @@ -1326,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. @@ -1595,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 @@ -1927,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)) @@ -2036,6 +2183,7 @@ end select continue close (isidep) + !----------------------------------------------------------------------- ! Calculate the "working" parameters of SC interactions. @@ -2144,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 @@ -2216,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 ! @@ -2310,11 +2551,11 @@ ! 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 - use control, only: rescode + use control, only: rescode,sugarcode ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.LOCAL' @@ -2327,7 +2568,7 @@ ! include 'COMMON.CONTROL' ! include 'COMMON.DISTFIT' ! include 'COMMON.SETUP' - integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,& + integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,& ! ishift_pdb logical :: lprn=.true.,fail real(kind=8),dimension(3) :: e1,e2,e3 @@ -2336,17 +2577,24 @@ character(len=5) :: atom character(len=80) :: card real(kind=8),dimension(3,20) :: sccor - integer :: kkk,lll,icha,kupa !rescode, + integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode, + integer :: isugar,molecprev + character*1 :: sugar real(kind=8) :: cou + real(kind=8),dimension(3) :: ccc !el local varables integer,dimension(2,maxres/3) :: hfrag_alloc integer,dimension(4,maxres/3) :: bfrag_alloc 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 ishift=0 + molecule=1 + counter=0 ! write (2,*) "UNRES_PDB",unres_pdb ires=0 ires_old=0 @@ -2384,9 +2632,11 @@ else if (card(:3).eq.'TER') then ! End current chain ires_old=ires+2 + ishift=ishift+1 ishift1=ishift1+1 - itype(ires_old,1)=ntyp1 - itype(ires_old-1,1)=ntyp1 + itype(ires_old,molecule)=ntyp1_molec(molecule) + itype(ires_old-1,molecule)=ntyp1_molec(molecule) + nres_molec(molecule)=nres_molec(molecule)+2 ibeg=2 ! write (iout,*) "Chain ended",ires,ishift,ires_old if (unres_pdb) then @@ -2402,6 +2652,7 @@ ! Read free energy if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp ! Fish out the ATOM cards. +! write(iout,*) 'card',card(1:20) if (index(card(1:4),'ATOM').gt.0) then read (card(12:16),*) atom ! write (iout,*) "! ",atom," !",ires @@ -2436,6 +2687,7 @@ if (res.ne.'GLY' .and. res.ne. 'ACE') then ishift=ishift-1 itype(1,1)=ntyp1 + nres_molec(molecule)=nres_molec(molecule)+1 endif ires=ires-ishift+ishift1 ires_old=ires @@ -2448,17 +2700,36 @@ ishift1=ishift1-1 !!!!! ! write (iout,*) "New chain started",ires,ishift,ishift1,"!" ires=ires-ishift+ishift1 +! print *,ires,ishift,ishift1 ires_old=ires ibeg=0 else ishift=ishift-(ires-ishift+ishift1-ires_old-1) ires=ires-ishift+ishift1 ires_old=ires - endif + endif +! print *,'atom',ires,atom if (res.eq.'ACE' .or. res.eq.'NHE') then itype(ires,1)=10 else - itype(ires,1)=rescode(ires,res,0,1) + if (atom.eq.'CA '.or.atom.eq.'N ') then + molecule=1 + itype(ires,molecule)=rescode(ires,res,0,molecule) +! nres_molec(molecule)=nres_molec(molecule)+1 + else + 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 ires=ires-ishift+ishift1 @@ -2471,31 +2742,101 @@ if (atom.eq.'CA' .or. atom.eq.'CH3' .or. & res.eq.'NHE'.and.atom(:2).eq.'HN') then read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) +! print *,ires,ishift,ishift1 ! write (iout,*) "backbone ",atom #ifdef DEBUG write (iout,'(2i3,2x,a,3f8.3)') & ires,itype(ires,1),res,(c(j,ires),j=1,3) #endif iii=iii+1 + nres_molec(molecule)=nres_molec(molecule)+1 do j=1,3 sccor(j,iii)=c(j,ires) enddo + else if (.not.unres_pdb .and. (atom.eq."C1'" .or. & + atom.eq."C2'" .or. atom.eq."C3'" & + .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 + counter=counter+1 +! iii=iii+1 +! do j=1,3 +! sccor(j,iii)=c(j,ires) +! enddo + 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 +! do j=1,3 +! sccor(j,iii)=c(j,ires) +! enddo + counter=0 + endif +! print *, "ATOM",atom(1:3) + 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) + else + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif ! write (*,*) card(23:27),ires,itype(ires,1) else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. & atom.ne.'N' .and. atom.ne.'C' .and. & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. & - atom.ne.'OXT' .and. atom(:2).ne.'3H') then + atom.ne.'OXT' .and. atom(:2).ne.'3H' & + .and. atom.ne.'P '.and. & + atom(1:1).ne.'H' .and. & + atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'& + ) then ! write (iout,*) "sidechain ",atom +! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3) + if ((molecule.ne.2).or.(atom(3:3).ne."'")) then +! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3) iii=iii+1 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif endif - endif + else if ((ions).and.(card(1:6).eq.'HETATM')) then + + read (card(12:16),*) atom + print *,"HETATOM", atom + read (card(18:20),'(a3)') res + if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.& + (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') & + .or.(atom(1:2).eq.'K ')) & + then + ires=ires+1 + if (molecule.ne.5) molecprev=molecule + molecule=5 + nres_molec(molecule)=nres_molec(molecule)+1 + itype(ires,molecule)=rescode(ires,res(2:4),0,molecule) + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) + endif + endif !atom enddo 10 write (iout,'(a,i5)') ' Number of residues found: ',ires if (ires.eq.0) return ! Calculate dummy residue coordinates inside the "chain" of a multichain ! system nres=ires + if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2 +! print *,'I have', nres_molec(:) + + do k=1,4 ! ions are without dummy + if (nres_molec(k).eq.0) cycle do i=2,nres-1 ! write (iout,*) i,itype(i,1) ! if (itype(i,1).eq.ntyp1) then @@ -2506,8 +2847,21 @@ ! dc(j,i)=c(j,i) ! enddo ! endif - if (itype(i,1).eq.ntyp1) then - if (itype(i+1,1).eq.ntyp1) then + 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" + do j=1,5 + if (itype(i+2,j).ne.0) then + itype(i+1,k)=0 + itype(i+1,j)=ntyp1_molec(j) + nres_molec(k)=nres_molec(k)-1 + nres_molec(j)=nres_molec(j)+1 + go to 3331 + endif + enddo + 3331 continue + endif ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false @@ -2520,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 @@ -2556,6 +2910,7 @@ endif !itype.eq.ntyp1 enddo + enddo ! Calculate the CM of the last side chain. if (iii.gt.0) then if (unres_pdb) then @@ -2569,9 +2924,12 @@ ! nres=ires nsup=nres nstart_sup=1 - if (itype(nres,1).ne.10) then +! print *,"molecule",molecule + if ((itype(nres,1).ne.10)) then nres=nres+1 - itype(nres,1)=ntyp1 + if (molecule.eq.5) molecule=molecprev + itype(nres,molecule)=ntyp1_molec(molecule) + nres_molec(molecule)=nres_molec(molecule)+1 if (unres_pdb) then ! 2/15/2013 by Adam: corrected insertion of the last dummy residue call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) @@ -2591,6 +2949,8 @@ enddo endif endif +! print *,'I have',nres, nres_molec(:) + !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb #ifdef WHAM_RUN if (nres.ne.nres0) then @@ -2638,6 +2998,7 @@ c(j,nres+1)=c(j,1) c(j,2*nres)=c(j,nres) enddo + if (itype(1,1).eq.ntyp1) then nsup=nsup-1 nstart_sup=2 @@ -2660,6 +3021,24 @@ enddo endif endif +! First lets assign correct dummy to correct type of chain +! 1) First residue + if (itype(1,1).eq.ntyp1) then + if (itype(2,1).eq.0) then + do j=2,5 + if (itype(2,j).ne.0) then + itype(1,1)=0 + itype(1,j)=ntyp1_molec(j) + nres_molec(1)=nres_molec(1)-1 + nres_molec(j)=nres_molec(j)+1 + go to 3231 + endif + enddo +3231 continue + endif + endif + print *,'I have',nres, nres_molec(:) + ! Copy the coordinates to reference coordinates ! do i=1,2*nres ! do j=1,3 @@ -2673,8 +3052,8 @@ write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" do ires=1,nres - write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') & - restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),& + write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') & + (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),& (c(j,ires+nres),j=1,3) enddo endif @@ -2689,7 +3068,75 @@ (c(j,nres+ires),j=1,3) enddo endif +! NOW LETS ROCK! SORTING + 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 + do i=1,nres + if (itype(i,k).ne.0) then + do j=1,3 + c_temporary(j,seqalingbegin)=c(j,i) + c_temporary(j,seqalingbegin+nres)=c(j,i+nres) + + enddo + itype_temporary(seqalingbegin,k)=itype(i,k) + istype_temp(seqalingbegin)=istype(i) + molnum(seqalingbegin)=k + seqalingbegin=seqalingbegin+1 + endif + enddo + enddo + do i=1,2*nres + do j=1,3 + c(j,i)=c_temporary(j,i) + enddo + enddo + 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 + nsup=nsup-1 + nstart_sup=2 + if (unres_pdb) then +! 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(2,3,4,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,1)=c(j,2)-1.9d0*e2(j) + enddo + else + do j=1,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 + endif + endif + if (lprn) then + write (iout,'(/a)') & + "Cartesian coordinates of the reference structure after sorting" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & + "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') & + (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),& + (c(j,ires+nres),j=1,3) + enddo + endif + +! print *,seqalingbegin,nres if(.not.allocated(vbld)) then allocate(vbld(2*nres)) do i=1,2*nres @@ -2768,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) @@ -3260,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.) @@ -3373,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+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)+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),& - gamsc(i),stdfsc(i)*dsqrt(gamsc(i)) + write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),& + gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1)) enddo endif else if (tbf) then @@ -3748,6 +4202,8 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',readonly,shared) + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared) call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',readonly,shared) call getenv_loc('ROTPAR',rotname) @@ -3762,6 +4218,19 @@ open (ielep,file=elename,status='old',readonly,shared) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly,shared) + + call getenv_loc('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared) + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared) + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old',readonly,shared) + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared) + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared) + + #elif (defined CRAY) || (defined AIX) open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& action='read') @@ -3772,6 +4241,9 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',action='read') + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old',action='read') + ! print *,"Processor",myrank," opened file IBOND" call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',action='read') @@ -3796,6 +4268,18 @@ ! print *,"Processor",myrank," opened file IELEP" call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',action='read') + + call getenv_loc('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old',action='read') + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old',action='read') + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old',action='read') + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old',action='read') + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old',action='read') + call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') call getenv_loc('TUBEPAR',tubename) @@ -3810,6 +4294,9 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old') + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old') + call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old') call getenv_loc('ROTPAR',rotname) @@ -3826,6 +4313,17 @@ open (ielep,file=elename,status='old') call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old') + + open (ithep_nucl,file=thetname_nucl,status='old') + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old') + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old') + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old') + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old') + call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old') call getenv_loc('TUBEPAR',tubename) @@ -3838,6 +4336,8 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',action='read') + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old',action='read') call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',action='read') call getenv_loc('ROTPAR',rotname) @@ -3860,6 +4360,18 @@ open (ielep,file=elename,status='old',readonly) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly) + + call getenv_loc('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old',action='read') + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old',action='read') + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old',action='read') + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old',action='read') + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old',action='read') + call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') call getenv_loc('TUBEPAR',tubename) @@ -3870,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