X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_config.f90;h=3f4303ac0a17a5e96f77b87baf4ce8af642b1d10;hb=705644e0cbb7678faefd6fe1bc436159d38ad85d;hp=0b10e1110166a63859d05a0e4ce433b459c6090e;hpb=35f220f409bd5d21be33a402d79da2c23d3e0c3a;p=unres4.git diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 0b10e11..3f4303a 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -5,6 +5,7 @@ use io_base use geometry_data use geometry + use control_data, only:maxterm_sccor implicit none !----------------------------------------------------------------------------- ! Max. number of residue types and parameters in expressions for @@ -21,7 +22,7 @@ ! parameter (maxtor=4,maxterm=10) !----------------------------------------------------------------------------- ! Max number of torsional terms in SCCOR - integer,parameter :: maxterm_sccor=6 +!el integer,parameter :: maxterm_sccor=6 !----------------------------------------------------------------------------- character(len=1),dimension(:),allocatable :: secstruc !(maxres) !----------------------------------------------------------------------------- @@ -29,7 +30,7 @@ ! !----------------------------------------------------------------------------- contains -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- ! bank.F io_csa !----------------------------------------------------------------------------- @@ -477,7 +478,7 @@ ! write (iout,100) ! do i=1,nres -! write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),& +! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),& ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i) ! enddo ! 100 format (//' alpha-carbon coordinates ',& @@ -575,15 +576,15 @@ endif enddo ndih_constr=ii - deallocate(secstruc) +! deallocate(secstruc) return 100 continue write(iout,'(A30,A80)')'Error reading file SECPRED',secpred - deallocate(secstruc) +! deallocate(secstruc) return 110 continue write(iout,'(A20)')'Error reading FTORS' - deallocate(secstruc) +! deallocate(secstruc) return end subroutine secstrp2dihc !----------------------------------------------------------------------------- @@ -689,7 +690,7 @@ ' in position',i4) return end subroutine read_secstr_pred -#endif +!#endif !----------------------------------------------------------------------------- ! parmread.F !----------------------------------------------------------------------------- @@ -697,7 +698,7 @@ use geometry_data use energy_data - use control_data, only:maxtor,maxterm + use control_data, only:maxterm !,maxtor use MD_data use MPI_data !el use map_data @@ -734,20 +735,20 @@ 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) :: buse character(len=3) :: lancuch !,ucase !el local variables integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm - integer :: maxinter,junk,kk,ii + integer :: maxinter,junk,kk,ii,ncatprotparm 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)) !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)) - ! ! For printing parameters after they are read set the following in the UNRES ! C-shell script: @@ -775,15 +776,23 @@ ! 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(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp) + allocate(long_r_sidechain(ntyp)) + allocate(short_r_sidechain(ntyp)) + dsc(:)=0.0d0 + dsc_inv(:)=0.0d0 + +#ifdef CRYST_BOND allocate(msc(ntyp+1)) !(ntyp+1) allocate(isc(ntyp+1)) !(ntyp+1) allocate(restok(ntyp+1)) !(ntyp+1) - allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp) -#ifdef CRYST_BOND read (ibond,*) vbldp0,akp,mp,ip,pstok do i=1,ntyp nbondterm(i)=1 @@ -796,10 +805,14 @@ endif enddo #else - read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok - do i=1,ntyp + allocate(msc(ntyp+1,5)) !(ntyp+1) + allocate(isc(ntyp+1,5)) !(ntyp+1) + allocate(restok(ntyp+1,5)) !(ntyp+1) + + 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 @@ -808,20 +821,57 @@ 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),nbondterm(i),& - vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i) + write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(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) enddo enddo endif + do i=1,ntyp_molec(5) + read(iion,*) msc(i,5),restok(i,5) + print *,msc(i,5),restok(i,5) + enddo + ip(5)=0.2 +! isc(5)=0.2 + read (iion,*) ncatprotparm + allocate(catprm(ncatprotparm,4)) + do k=1,4 + read (iion,*) (catprm(i,k),i=1,ncatprotparm) + enddo + print *, catprm +! read (iion,*) (vcatprm(k),k=1,ncatprotpram) !---------------------------------------------------- allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp)) allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp) @@ -829,26 +879,26 @@ allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1) allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp) allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp) - do i=-ntyp,ntyp - a0thet(i)=0.0D0 - do j=1,2 - do ichir1=-1,1 - do ichir2=-1,1 - athet(j,i,ichir1,ichir2)=0.0D0 - bthet(j,i,ichir1,ichir2)=0.0D0 - enddo - enddo - enddo - do j=0,3 - polthet(j,i)=0.0D0 - enddo - do j=1,3 - gthet(j,i)=0.0D0 - enddo - theta0(i)=0.0D0 - sig0(i)=0.0D0 - sigc0(i)=0.0D0 - enddo + + a0thet(:)=0.0D0 + athet(:,:,:,:)=0.0D0 + bthet(:,:,:,:)=0.0D0 + polthet(:,:)=0.0D0 + gthet(:,:)=0.0D0 + 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 ! @@ -907,7 +957,7 @@ ' ATHETA0 ',' A1 ',' A2 ',& ' B1 ',' B2 ' do i=1,ntyp - write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,& + write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,& a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2) enddo write (iout,'(/a/9x,5a/79(1h-))') & @@ -915,7 +965,7 @@ ' ALPH0 ',' ALPH1 ',' ALPH2 ',& ' ALPH3 ',' SIGMA0C ' do i=1,ntyp - write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,& + write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,& (polthet(j,i),j=0,3),sigc0(i) enddo write (iout,'(/a/9x,5a/79(1h-))') & @@ -923,7 +973,7 @@ ' THETA0 ',' SIGMA0 ',' G1 ',& ' G2 ',' G3 ' do i=1,ntyp - write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),& + write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),& sig0(i),(gthet(j,i),j=1,3) enddo else @@ -934,7 +984,7 @@ ' theta0 ',' a1*10^2 ',' a2*10^2 ',& ' b1*10^1 ',' b2*10^1 ' do i=1,ntyp - write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),& + write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),& a0thet(i),(100*athet(j,i,1,1),j=1,2),& (10*bthet(j,i,1,1),j=1,2) enddo @@ -943,7 +993,7 @@ ' alpha0 ',' alph1 ',' alph2 ',& ' alhp3 ',' sigma0c ' do i=1,ntyp - write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),& + write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),& (polthet(j,i),j=0,3),sigc0(i) enddo write (iout,'(/a/9x,5a/79(1h-))') & @@ -951,7 +1001,7 @@ ' theta0 ',' sigma0*10^2 ',' G1*10^-1',& ' G2 ',' G3*10^1 ' do i=1,ntyp - write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),& + write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),& 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0 enddo endif @@ -967,27 +1017,27 @@ !---------------------------------------------------- allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) - allocate(aa0thet(-maxthetyp1:maxthetyp1,& - -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) + allocate(aa0thet(-nthetyp-1:nthetyp+1,& + -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) - allocate(aathet(ntheterm,-maxthetyp1:maxthetyp1,& - -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) + allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,& + -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(maxtheterm,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) - allocate(bbthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,& - -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) - allocate(ccthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,& - -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) - allocate(ddthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,& - -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) - allocate(eethet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,& - -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) + allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& + -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) + allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& + -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) + allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& + -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) + allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& + -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) - allocate(ffthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,& - -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) - allocate(ggthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,& - -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) + allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,& + -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) + allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,& + -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) @@ -995,34 +1045,16 @@ do i=-ntyp1,-1 ithetyp(i)=-ithetyp(-i) enddo - do iblock=1,2 - do i=-maxthetyp,maxthetyp - do j=-maxthetyp,maxthetyp - do k=-maxthetyp,maxthetyp - aa0thet(i,j,k,iblock)=0.0d0 - do l=1,ntheterm - aathet(l,i,j,k,iblock)=0.0d0 - enddo - do l=1,ntheterm2 - do m=1,nsingle - bbthet(m,l,i,j,k,iblock)=0.0d0 - ccthet(m,l,i,j,k,iblock)=0.0d0 - ddthet(m,l,i,j,k,iblock)=0.0d0 - eethet(m,l,i,j,k,iblock)=0.0d0 - enddo - enddo - do l=1,ntheterm3 - do m=1,ndouble - do mm=1,ndouble - ffthet(mm,m,l,i,j,k,iblock)=0.0d0 - ggthet(mm,m,l,i,j,k,iblock)=0.0d0 - enddo - enddo - enddo - enddo - enddo - enddo - enddo + + aa0thet(:,:,:,:)=0.0d0 + aathet(:,:,:,:,:)=0.0d0 + bbthet(:,:,:,:,:,:)=0.0d0 + ccthet(:,:,:,:,:,:)=0.0d0 + ddthet(:,:,:,:,:,:)=0.0d0 + eethet(:,:,:,:,:,:)=0.0d0 + ffthet(:,:,:,:,:,:,:)=0.0d0 + ggthet(:,:,:,:,:,:,:)=0.0d0 + ! VAR:iblock means terminally blocking group 1=non-proline 2=proline do iblock=1,2 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine @@ -1220,6 +1252,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) @@ -1227,28 +1318,13 @@ allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp) allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp) - do i=1,ntyp - do j=1,maxlob - bsc(j,i)=0.0D0 - nlob(i)=0 - enddo - enddo - nlob(ntyp1)=0 - dsc(ntyp1)=0.0D0 - - do i=-ntyp,ntyp - do j=1,maxlob - do k=1,3 - censc(k,j,i)=0.0D0 - enddo - do k=1,3 - do l=1,3 - gaussc(l,k,j,i)=0.0D0 - enddo - enddo - enddo - enddo - + bsc(:,:)=0.0D0 + nlob(:)=0 + nlob(:)=0 + dsc(:)=0.0D0 + censc(:,:,:)=0.0D0 + gaussc(:,:,:,:)=0.0D0 + #ifdef CRYST_SC ! ! Read the parameters of the probability distribution/energy expression @@ -1313,7 +1389,7 @@ nlobi=nlob(i) if (nlobi.gt.0) then if (LaTeX) then - write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),& + write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),& ' # of gaussian lobes:',nlobi,' dsc:',dsc(i) write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') & 'log h',(bsc(j,i),j=1,nlobi) @@ -1354,6 +1430,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. @@ -1451,22 +1545,18 @@ allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) !el--------------------------- - do iblock=1,2 - do i=-ntortyp,ntortyp - do j=-ntortyp,ntortyp - nterm(i,j,iblock)=0 - nlor(i,j,iblock)=0 - enddo - enddo - enddo + nterm(:,:,:)=0 + nlor(:,:,:)=0 !el--------------------------- read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp) do i=-ntyp,-1 itortyp(i)=-itortyp(-i) enddo - write (iout,*) 'ntortyp',ntortyp + itortyp(ntyp1)=ntortyp + itortyp(-ntyp1)=-ntortyp do iblock=1,2 + write (iout,*) 'ntortyp',ntortyp do i=0,ntortyp-1 do j=-ntortyp+1,ntortyp-1 read (itorp,*,end=113,err=113) nterm(i,j,iblock),& @@ -1543,7 +1633,7 @@ if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) & .or. t3.ne.toronelet(k)) then - write (iout,*) "Error in double torsional parameter file",& + write (iout,*) "Error in double torsional parameter file",& i,j,k,t1,t2,t3 #ifdef MPI call MPI_Finalize(Ierror) @@ -1627,6 +1717,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 @@ -1770,6 +1902,7 @@ #endif if (lprint) then + l=3 write (iout,'(/a/)') 'Torsional constants:' do i=1,nsccortyp do j=1,nsccortyp @@ -1797,7 +1930,7 @@ write (iout,*) "Coefficients of the cumulants" endif read (ifourier,*) nloctyp -write(iout,*) "nloctyp",nloctyp +!write(iout,*) "nloctyp",nloctyp !el from module energy------- allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor) allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor) @@ -1807,92 +1940,95 @@ write(iout,*) "nloctyp",nloctyp allocate(ee(2,2,-nloctyp-1:nloctyp+1)) allocate(ctilde(2,2,-nloctyp-1:nloctyp+1)) allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor) +! el + b1(1,:)=0.0d0 + b1(2,:)=0.0d0 !-------------------------------- 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) (buse(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,*) 'Type',i + write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(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)=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)=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)=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) + B1(1,i) = buse(3) + B1(2,i) = buse(5) + B1(1,-i) = buse(3) + B1(2,-i) = -buse(5) +! buse1(1,i)=0.0d0 +! buse1(2,i)=0.0d0 + B1tilde(1,i) = buse(3) + B1tilde(2,i) =-buse(5) + B1tilde(1,-i) =-buse(3) + B1tilde(2,-i) =buse(5) +! buse1tilde(1,i)=0.0d0 +! buse1tilde(2,i)=0.0d0 + B2(1,i) = buse(2) + B2(2,i) = buse(4) + B2(1,-i) =buse(2) + B2(2,-i) =-buse(4) + +! buse2(1,i)=0.0d0 +! buse2(2,i)=0.0d0 + CC(1,1,i)= buse(7) + CC(2,2,i)=-buse(7) + CC(2,1,i)= buse(9) + CC(1,2,i)= buse(9) + CC(1,1,-i)= buse(7) + CC(2,2,-i)=-buse(7) + CC(2,1,-i)=-buse(9) + CC(1,2,-i)=-buse(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)=buse(7) + Ctilde(1,2,i)=buse(9) + Ctilde(2,1,i)=-buse(9) + Ctilde(2,2,i)=buse(7) + Ctilde(1,1,-i)=buse(7) + Ctilde(1,2,-i)=-buse(9) + Ctilde(2,1,-i)=buse(9) + Ctilde(2,2,-i)=buse(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)= buse(6) + DD(2,2,i)=-buse(6) + DD(2,1,i)= buse(8) + DD(1,2,i)= buse(8) + DD(1,1,-i)= buse(6) + DD(2,2,-i)=-buse(6) + DD(2,1,-i)=-buse(8) + DD(1,2,-i)=-buse(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)=buse(6) + Dtilde(1,2,i)=buse(8) + Dtilde(2,1,i)=-buse(8) + Dtilde(2,2,i)=buse(6) + Dtilde(1,1,-i)=buse(6) + Dtilde(1,2,-i)=-buse(8) + Dtilde(2,1,-i)=buse(8) + Dtilde(2,2,-i)=buse(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)= buse(10)+buse(11) + EE(2,2,i)=-buse(10)+buse(11) + EE(2,1,i)= buse(12)-buse(13) + EE(1,2,i)= buse(12)+buse(13) + EE(1,1,-i)= buse(10)+buse(11) + EE(2,2,-i)=-buse(10)+buse(11) + EE(2,1,-i)=-buse(12)+buse(13) + EE(1,2,-i)=-buse(12)-buse(13) ! ee(1,1,i)=1.0d0 ! ee(2,2,i)=1.0d0 @@ -1921,7 +2057,6 @@ write(iout,*) "nloctyp",nloctyp enddo enddo endif - ! ! Read electrostatic-interaction parameters ! @@ -1957,18 +2092,17 @@ write(iout,*) "nloctyp",nloctyp 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) - do i=1,ntyp - do j=1,ntyp - augm(i,j)=0.0D0 - enddo - chip(i)=0.0D0 - alp(i)=0.0D0 - sigma0(i)=0.0D0 - sigii(i)=0.0D0 - rr0(i)=0.0D0 - enddo + allocate(epslip(ntyp,ntyp)) + augm(:,:)=0.0D0 + chip(:)=0.0D0 + alp(:)=0.0D0 + sigma0(:)=0.0D0 + sigii(:)=0.0D0 + rr0(:)=0.0D0 + !-------------------------------- read (isidep,*,end=117,err=117) ipot,expon @@ -1984,93 +2118,122 @@ write(iout,*) "nloctyp",nloctyp if(me.eq.king) & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),& ', exponents are ',expon,2*expon - goto (10,20,30,30,40) ipot +! goto (10,20,30,30,40) ipot + select case(ipot) !----------------------- LJ potential --------------------------------- - 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),& - (sigma0(i),i=1,ntyp) - if (lprint) then - write (iout,'(/a/)') 'Parameters of the LJ potential:' - write (iout,'(a/)') 'The epsilon array:' - call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) - write (iout,'(/a)') 'One-body parameters:' - write (iout,'(a,4x,a)') 'residue','sigma' - write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp) - endif - goto 50 + case (1) +! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),& + read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),& + (sigma0(i),i=1,ntyp) + if (lprint) then + write (iout,'(/a/)') 'Parameters of the LJ potential:' + write (iout,'(a/)') 'The epsilon array:' + call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) + write (iout,'(/a)') 'One-body parameters:' + write (iout,'(a,4x,a)') 'residue','sigma' + write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp) + endif +! goto 50 !----------------------- LJK potential -------------------------------- - 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),& - (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp) - if (lprint) then - write (iout,'(/a/)') 'Parameters of the LJK potential:' - write (iout,'(a/)') 'The epsilon array:' - call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) - write (iout,'(/a)') 'One-body parameters:' - write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 ' - write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),& - i=1,ntyp) - endif - goto 50 + case(2) +! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),& + read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),& + (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp) + if (lprint) then + write (iout,'(/a/)') 'Parameters of the LJK potential:' + write (iout,'(a/)') 'The epsilon array:' + call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) + write (iout,'(/a)') 'One-body parameters:' + write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 ' + write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),& + i=1,ntyp) + endif +! goto 50 !---------------------- GB or BP potential ----------------------------- - 30 do i=1,ntyp - read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp) - enddo - read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp) - 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) -! For the GB potential convert sigma'**2 into chi' - if (ipot.eq.4) then - do i=1,ntyp - chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0) + case(3:4) +! 30 do i=1,ntyp + do i=1,ntyp + read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp) enddo - endif - if (lprint) then - write (iout,'(/a/)') 'Parameters of the BP potential:' - write (iout,'(a/)') 'The epsilon array:' - call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) - write (iout,'(/a)') 'One-body parameters:' - write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',& - ' chip ',' alph ' - write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),& - chip(i),alp(i),i=1,ntyp) - endif - goto 50 + read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp) + 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 + chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0) + enddo + endif + if (lprint) then + write (iout,'(/a/)') 'Parameters of the BP potential:' + write (iout,'(a/)') 'The epsilon array:' + call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) + write (iout,'(/a)') 'One-body parameters:' + write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',& + ' chip ',' alph ' + write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),& + chip(i),alp(i),i=1,ntyp) + endif +! goto 50 !--------------------- GBV potential ----------------------------------- - 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),& - (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),& - (chip(i),i=1,ntyp),(alp(i),i=1,ntyp) - if (lprint) then - write (iout,'(/a/)') 'Parameters of the GBV potential:' - write (iout,'(a/)') 'The epsilon array:' - call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) - write (iout,'(/a)') 'One-body parameters:' - write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',& - 's||/s_|_^2',' chip ',' alph ' - write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),& - sigii(i),chip(i),alp(i),i=1,ntyp) - endif - 50 continue + case(5) +! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),& + read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),& + (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),& + (chip(i),i=1,ntyp),(alp(i),i=1,ntyp) + if (lprint) then + write (iout,'(/a/)') 'Parameters of the GBV potential:' + write (iout,'(a/)') 'The epsilon array:' + call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) + write (iout,'(/a)') 'One-body parameters:' + write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',& + 's||/s_|_^2',' chip ',' alph ' + write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),& + sigii(i),chip(i),alp(i),i=1,ntyp) + endif + case default + write(iout,*)"Wrong ipot" +! 50 continue + end select + continue close (isidep) + !----------------------------------------------------------------------- ! 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) - do i=1,ntyp1 - do j=1,ntyp1 - aa(i,j)=0.0D0 - bb(i,j)=0.0D0 - chi(i,j)=0.0D0 - sigma(i,j)=0.0D0 - r0(i,j)=0.0D0 - enddo - enddo + 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) + eps(i,j)=eps(j,i) + epslip(i,j)=epslip(j,i) enddo enddo do i=1,ntyp @@ -2099,10 +2262,18 @@ write(iout,*) "nloctyp",nloctyp 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 @@ -2135,18 +2306,280 @@ write(iout,*) "nloctyp",nloctyp 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,1),restyp(j,1),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 - allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2) - do i=1,ntyp - do j=1,2 - bad(i,j)=0.0D0 + 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 + 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 +!-----------------READING SC BASE POTENTIALS----------------------------- + allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2)) + allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2)) + allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2))) + allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2)) + allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2))) + allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2))) + allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2))) + allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2))) + allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2))) + allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2))) + allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2))) + + + do i=1,ntyp_molec(1) + do j=1,ntyp_molec(2)-1 ! without U then we will take T for U + write (*,*) "Im in ", i, " ", j + read(isidep_scbase,*) & + eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),& + chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2) + write(*,*) "eps",eps_scbase(i,j) + read(isidep_scbase,*) & + (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), & + chis_scbase(i,j,1),chis_scbase(i,j,2) + read(isidep_scbase,*) & + dhead_scbasei(i,j), & + dhead_scbasej(i,j), & + rborn_scbasei(i,j),rborn_scbasej(i,j) + read(isidep_scbase,*) & + (wdipdip_scbase(k,i,j),k=1,3), & + (wqdip_scbase(k,i,j),k=1,2) + read(isidep_scbase,*) & + alphapol_scbase(i,j), & + epsintab_scbase(i,j) + END DO + END DO + allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2))) + + do i=1,ntyp_molec(1) + do j=1,ntyp_molec(2)-1 + epsij=eps_scbase(i,j) + rrij=sigma_scbase(i,j) +! r0(i,j)=rrij +! r0(j,i)=rrij + rrij=rrij**expon +! epsij=eps(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_scbase(i,j)=epsij*rrij*rrij + bb_scbase(i,j)=-sigeps*epsij*rrij + enddo + enddo +!-----------------READING PEP BASE POTENTIALS------------------- + allocate(eps_pepbase(ntyp_molec(2))) + allocate(sigma_pepbase(ntyp_molec(2))) + allocate(chi_pepbase(ntyp_molec(2),2)) + allocate(chipp_pepbase(ntyp_molec(2),2)) + allocate(alphasur_pepbase(4,ntyp_molec(2))) + allocate(sigmap1_pepbase(ntyp_molec(2))) + allocate(sigmap2_pepbase(ntyp_molec(2))) + allocate(chis_pepbase(ntyp_molec(2),2)) + allocate(wdipdip_pepbase(3,ntyp_molec(2))) + + + do j=1,ntyp_molec(2)-1 ! without U then we will take T for U + write (*,*) "Im in ", i, " ", j + read(isidep_pepbase,*) & + eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),& + chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2) + write(*,*) "eps",eps_pepbase(j) + read(isidep_pepbase,*) & + (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), & + chis_pepbase(j,1),chis_pepbase(j,2) + read(isidep_pepbase,*) & + (wdipdip_pepbase(k,j),k=1,3) + END DO + allocate(aa_pepbase(ntyp_molec(2))) + allocate(bb_pepbase(ntyp_molec(2))) + + do j=1,ntyp_molec(2)-1 + epsij=eps_pepbase(j) + rrij=sigma_pepbase(j) +! r0(i,j)=rrij +! r0(j,i)=rrij + rrij=rrij**expon +! epsij=eps(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_pepbase(j)=epsij*rrij*rrij + bb_pepbase(j)=-sigeps*epsij*rrij + enddo +!--------------READING SC PHOSPHATE------------------------------------- + allocate(eps_scpho(ntyp_molec(1))) + allocate(sigma_scpho(ntyp_molec(1))) + allocate(chi_scpho(ntyp_molec(1),2)) + allocate(chipp_scpho(ntyp_molec(1),2)) + allocate(alphasur_scpho(4,ntyp_molec(1))) + allocate(sigmap1_scpho(ntyp_molec(1))) + allocate(sigmap2_scpho(ntyp_molec(1))) + allocate(chis_scpho(ntyp_molec(1),2)) + allocate(wqq_scpho(ntyp_molec(1))) + allocate(wqdip_scpho(2,ntyp_molec(1))) + allocate(alphapol_scpho(ntyp_molec(1))) + allocate(epsintab_scpho(ntyp_molec(1))) + allocate(dhead_scphoi(ntyp_molec(1))) + allocate(rborn_scphoi(ntyp_molec(1))) + allocate(rborn_scphoj(ntyp_molec(1))) + allocate(alphi_scpho(ntyp_molec(1))) + + +! j=1 + do j=1,ntyp_molec(1) ! without U then we will take T for U + write (*,*) "Im in scpho ", i, " ", j + read(isidep_scpho,*) & + eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),& + chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2) + write(*,*) "eps",eps_scpho(j) + read(isidep_scpho,*) & + (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), & + chis_scpho(j,1),chis_scpho(j,2) + read(isidep_scpho,*) & + (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j) + read(isidep_scpho,*) & + epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), & + alphi_scpho(j) + + END DO + allocate(aa_scpho(ntyp_molec(1))) + allocate(bb_scpho(ntyp_molec(1))) + + do j=1,ntyp_molec(1) + epsij=eps_scpho(j) + rrij=sigma_scpho(j) +! r0(i,j)=rrij +! r0(j,i)=rrij + rrij=rrij**expon +! epsij=eps(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_scpho(j)=epsij*rrij*rrij + bb_scpho(j)=-sigeps*epsij*rrij + enddo + + + read(isidep_peppho,*) & + eps_peppho,sigma_peppho + read(isidep_peppho,*) & + (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho + read(isidep_peppho,*) & + (wqdip_peppho(k),k=1,2) + + epsij=eps_peppho + rrij=sigma_peppho +! r0(i,j)=rrij +! r0(j,i)=rrij + rrij=rrij**expon +! epsij=eps(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_peppho=epsij*rrij*rrij + bb_peppho=-sigeps*epsij*rrij + + + allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2) + bad(:,:)=0.0D0 #ifdef OLDSCP ! @@ -2197,6 +2630,20 @@ write(iout,*) "nloctyp",nloctyp 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 ! @@ -2259,88 +2706,7 @@ write(iout,*) "nloctyp",nloctyp stop return end subroutine parmread -!----------------------------------------------------------------------------- -! permut.F -!----------------------------------------------------------------------------- - subroutine permut(isym) - - use geometry_data, only: tabperm -! use energy_data -! use control_data, only:lsecondary -! use MD_data -! use MPI_data -! use map_data -! use energy -! use geometry -! use control -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.LOCAL' -! include 'COMMON.VAR' -! include 'COMMON.CHAIN' -! include 'COMMON.INTERACT' -! include 'COMMON.IOUNITS' -! include 'COMMON.GEO' -! include 'COMMON.NAMES' -! include 'COMMON.CONTROL' - - integer :: n,isym -! logical nextp -!el external nextp - integer,dimension(isym) :: a -! parameter(n=symetr) -!el local variables - integer :: kkk,i - - n=isym - if (n.eq.1) then - tabperm(1,1)=1 - return - endif - kkk=0 - do i=1,n - a(i)=i - enddo - 10 print *,(a(i),i=1,n) - kkk=kkk+1 - do i=1,n - tabperm(kkk,i)=a(i) -! write (iout,*) "tututu", kkk - enddo - if(nextp(n,a)) go to 10 - return - end subroutine permut -!----------------------------------------------------------------------------- - logical function nextp(n,a) - - integer :: n,i,j,k,t -! logical :: nextp - integer,dimension(n) :: a - i=n-1 - 10 if(a(i).lt.a(i+1)) go to 20 - i=i-1 - if(i.eq.0) go to 20 - go to 10 - 20 j=i+1 - k=n - 30 t=a(j) - a(j)=a(k) - a(k)=t - j=j+1 - k=k-1 - if(j.lt.k) go to 30 - j=i - if(j.ne.0) go to 40 - nextp=.false. - return - 40 j=j+1 - if(a(j).lt.a(i)) go to 40 - t=a(i) - a(i)=a(j) - a(j)=t - nextp=.true. - return - end function nextp +#endif !----------------------------------------------------------------------------- ! printmat.f !----------------------------------------------------------------------------- @@ -2372,11 +2738,11 @@ write(iout,*) "nloctyp",nloctyp ! 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' @@ -2389,7 +2755,7 @@ write(iout,*) "nloctyp",nloctyp ! 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 @@ -2398,58 +2764,78 @@ write(iout,*) "nloctyp",nloctyp 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,firstion + 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 +#ifdef WHAM_RUN + do i=1,nres + do j=1,5 + itype(i,j)=0 + enddo + enddo +#endif nres=0 iii=0 lsecondary=.false. nhfrag=0 nbfrag=0 + do j=1,5 + nres_molec(j)=0 + enddo + + !----------------------------- allocate(hfrag(2,maxres/3)) !(2,maxres/3) allocate(bfrag(4,maxres/3)) !(4,maxres/3) - -!elwrite(iout,*)"poczatek read pdb" - + if(.not. allocated(istype)) allocate(istype(maxres)) do i=1,100000 read (ipdbin,'(a80)',end=10) card ! write (iout,'(a)') card if (card(:5).eq.'HELIX') then - nhfrag=nhfrag+1 - lsecondary=.true. - read(card(22:25),*) hfrag(1,nhfrag) - read(card(34:37),*) hfrag(2,nhfrag) + nhfrag=nhfrag+1 + lsecondary=.true. + read(card(22:25),*) hfrag(1,nhfrag) + read(card(34:37),*) hfrag(2,nhfrag) endif if (card(:5).eq.'SHEET') then - nbfrag=nbfrag+1 - lsecondary=.true. - read(card(24:26),*) bfrag(1,nbfrag) - read(card(35:37),*) bfrag(2,nbfrag) + nbfrag=nbfrag+1 + lsecondary=.true. + read(card(24:26),*) bfrag(1,nbfrag) + read(card(35:37),*) bfrag(2,nbfrag) !rc---------------------------------------- !rc to be corrected !!! - bfrag(3,nbfrag)=bfrag(1,nbfrag) - bfrag(4,nbfrag)=bfrag(2,nbfrag) + bfrag(3,nbfrag)=bfrag(1,nbfrag) + bfrag(4,nbfrag)=bfrag(2,nbfrag) !rc---------------------------------------- endif if (card(:3).eq.'END') then goto 10 else if (card(:3).eq.'TER') then ! End current chain - ires_old=ires+1 + ires_old=ires+2 + ishift=ishift+1 ishift1=ishift1+1 - itype(ires_old)=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 @@ -2458,12 +2844,14 @@ write(iout,*) "nloctyp",nloctyp enddo else call sccenter(ires,iii,sccor) +! iii=0 endif iii=0 endif ! 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 @@ -2481,7 +2869,7 @@ write(iout,*) "nloctyp",nloctyp ! write (iout,*) "Calculating sidechain center iii",iii if (unres_pdb) then do j=1,3 - dc(j,ires+nres)=sccor(j,iii) + dc(j,ires+ishift1-ishift-1)=sccor(j,iii) enddo else call sccenter(ires_old,iii,sccor) @@ -2497,30 +2885,51 @@ write(iout,*) "nloctyp",nloctyp ishift=ires-1 if (res.ne.'GLY' .and. res.ne. 'ACE') then ishift=ishift-1 - itype(1)=ntyp1 + itype(1,1)=ntyp1 + nres_molec(molecule)=nres_molec(molecule)+1 endif ires=ires-ishift+ishift1 ires_old=ires -! write (iout,*) "ishift",ishift," ires",ires, -! & " ires_old",ires_old +! write (iout,*) "ishift",ishift," ires",ires,& +! " ires_old",ires_old ibeg=0 else if (ibeg.eq.2) then ! Start a new chain -! ishift=-ires_old+ires-1 -! ishift1=ishift1+1 + ishift=-ires_old+ires-1 !!!!! + 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)=10 + itype(ires,1)=10 + else + if (atom.eq.'CA '.or.atom.eq.'N ') then + molecule=1 + itype(ires,molecule)=rescode(ires,res,0,molecule) + firstion=0 +! nres_molec(molecule)=nres_molec(molecule)+1 else - itype(ires)=rescode(ires,res,0) + molecule=2 + itype(ires,molecule)=rescode(ires,res(2:3),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 @@ -2529,45 +2938,200 @@ write(iout,*) "nloctyp",nloctyp if (card(27:27).eq."A" .or. card(27:27).eq."B") then ! ishift1=ishift1+1 endif -! write (2,*) "ires",ires," res ",res," ity",ity +! write (2,*) "ires",ires," res ",res!," ity"!,ity 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) -! write (iout,*) "backbone ",atom +! print *,ires,ishift,ishift1 +! write (iout,*) "backbone ",atom #ifdef DEBUG write (iout,'(2i3,2x,a,3f8.3)') & - ires,itype(ires),res,(c(j,ires),j=1,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 -! write (*,*) card(23:27),ires,itype(ires) + 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 + firstion=0 +! 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 + molecule=2 +! print *,"nres_molec(molecule)",nres_molec(molecule),ires + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) + nres_molec(molecule)=nres_molec(molecule)+1 + print *,"nres_molec(molecule)",nres_molec(molecule),ires + + else + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif + else if ((atom.eq."C1'").and.unres_pdb) then + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) +! 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 + if (firstion.eq.0) then + firstion=1 + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + endif + 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 + print *,"HERE",nres_molec(molecule) + res=res(2:3)//' ' + itype(ires,molecule)=rescode(ires,res,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).and.(molecule.ne.5)) & + .and.(.not.unres_pdb)) & + nres_molec(molecule)=nres_molec(molecule)-2 + print *,'I have',nres, 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) - if (itype(i).eq.ntyp1) then -! write (iout,*) "dummy",i,itype(i) - do j=1,3 - c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2 +! write (iout,*) i,itype(i,1) +! if (itype(i,1).eq.ntyp1) then +! write (iout,*) "dummy",i,itype(i,1) +! do j=1,3 +! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2 - dc(j,i)=c(j,i) - enddo - endif +! dc(j,i)=c(j,i) +! enddo +! endif + 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 +! if (unres_pdb) then +! 2/15/2013 by Adam: corrected insertion of the last dummy residue +! print *,i,'tu dochodze' +! call refsys(i-3,i-2,i-1,e1,e2,e3,fail) +! if (fail) then +! e2(1)=0.0d0 +! e2(2)=1.0d0 +! e2(3)=0.0d0 +! endif !fail +! print *,i,'a tu?' +! do j=1,3 +! c(j,i)=c(j,i-1)-1.9d0*e2(j) +! enddo +! else !unres_pdb + do j=1,3 + dcj=(c(j,i-2)-c(j,i-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i-1)+dcj + c(j,nres+i)=c(j,i) + enddo +! endif !unres_pdb + else !itype(i+1,1).eq.ntyp1 +! if (unres_pdb) then +! 2/15/2013 by Adam: corrected insertion of the first dummy residue +! call refsys(i+1,i+2,i+3,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,i)=c(j,i+1)-1.9d0*e2(j) + enddo +! else !unres_pdb + do j=1,3 + dcj=(c(j,i+3)-c(j,i+2))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i+1)-dcj + c(j,nres+i)=c(j,i) + enddo +! endif !unres_pdb + endif !itype(i+1,1).eq.ntyp1 + endif !itype.eq.ntyp1 + + enddo enddo ! Calculate the CM of the last side chain. if (iii.gt.0) then @@ -2582,56 +3146,69 @@ write(iout,*) "nloctyp",nloctyp ! nres=ires nsup=nres nstart_sup=1 - if (itype(nres).ne.10) then +! print *,"molecule",molecule + if ((itype(nres,1).ne.10)) then nres=nres+1 - itype(nres)=ntyp1 - if (unres_pdb) then + 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) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif - do j=1,3 - c(j,nres)=c(j,nres-1)-3.8d0*e2(j) - enddo - else +! call refsys(nres-3,nres-2,nres-1,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,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 - endif +! 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 + write (iout,*) "Error: wrong parameter value: NRES=",nres,& + " NRES0=",nres0 + stop "Error nres value in WHAM input" + endif +#endif !--------------------------------- !el reallocate tables - do i=1,maxres/3 - do j=1,2 - hfrag_alloc(j,i)=hfrag(j,i) - enddo - do j=1,4 - bfrag_alloc(j,i)=bfrag(j,i) - enddo - enddo +! do i=1,maxres/3 +! do j=1,2 +! hfrag_alloc(j,i)=hfrag(j,i) +! enddo +! do j=1,4 +! bfrag_alloc(j,i)=bfrag(j,i) +! enddo +! enddo - deallocate(hfrag) - deallocate(bfrag) - allocate(hfrag(2,nres/3)) !(2,maxres/3) +! deallocate(hfrag) +! deallocate(bfrag) +! allocate(hfrag(2,nres/3)) !(2,maxres/3) !el allocate(hfrag(2,nhfrag)) !(2,maxres/3) !el allocate(bfrag(4,nbfrag)) !(4,maxres/3) - allocate(bfrag(4,nres/3)) !(4,maxres/3) +! allocate(bfrag(4,nres/3)) !(4,maxres/3) - do i=1,nhfrag - do j=1,2 - hfrag(j,i)=hfrag_alloc(j,i) - enddo - enddo - do i=1,nbfrag - do j=1,4 - bfrag(j,i)=bfrag_alloc(j,i) - enddo - enddo +! do i=1,nhfrag +! do j=1,2 +! hfrag(j,i)=hfrag_alloc(j,i) +! enddo +! enddo +! do i=1,nbfrag +! do j=1,4 +! bfrag(j,i)=bfrag_alloc(j,i) +! enddo +! enddo !el end reallocate tables !--------------------------------- do i=2,nres-1 @@ -2643,7 +3220,8 @@ write(iout,*) "nloctyp",nloctyp c(j,nres+1)=c(j,1) c(j,2*nres)=c(j,nres) enddo - if (itype(1).eq.ntyp1) then + + if (itype(1,1).eq.ntyp1) then nsup=nsup-1 nstart_sup=2 if (unres_pdb) then @@ -2655,16 +3233,34 @@ write(iout,*) "nloctyp",nloctyp 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 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 @@ -2678,29 +3274,115 @@ write(iout,*) "nloctyp",nloctyp 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)),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 +! znamy już nres wiec mozna alokowac tablice ! Calculate internal coordinates. if(me.eq.king.or..not.out1file)then write (iout,'(a)') & "Backbone and SC coordinates as read from the PDB" do ires=1,nres write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') & - ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),& + ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),& (c(j,nres+ires),j=1,3) enddo endif +! NOW LETS ROCK! SORTING + allocate(c_temporary(3,2*nres)) + allocate(itype_temporary(nres,5)) + if (.not.allocated(molnum)) allocate(molnum(nres+1)) + if (.not.allocated(istype)) write(iout,*) & + "SOMETHING WRONG WITH ISTYTPE" + 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) + print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin + 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 + vbld(i)=0.d0 + enddo + endif + if(.not.allocated(vbld_inv)) then + allocate(vbld_inv(2*nres)) + do i=1,2*nres + vbld_inv(i)=0.d0 + enddo + endif +!!!el + if(.not.allocated(theta)) then + allocate(theta(nres+2)) + theta(:)=0.0d0 + endif - if(.not.allocated(vbld)) allocate(vbld(2*nres)) - if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres)) - if(.not.allocated(theta)) allocate(theta(nres+2)) if(.not.allocated(phi)) allocate(phi(nres+2)) if(.not.allocated(alph)) allocate(alph(nres+2)) if(.not.allocated(omeg)) allocate(omeg(nres+2)) - if(.not.allocated(theta)) allocate(theta(nres+2)) if(.not.allocated(thetaref)) allocate(thetaref(nres+2)) if(.not.allocated(phiref)) allocate(phiref(nres+2)) if(.not.allocated(costtab)) allocate(costtab(nres)) @@ -2710,23 +3392,29 @@ write(iout,*) "nloctyp",nloctyp if(.not.allocated(xxref)) allocate(xxref(nres)) if(.not.allocated(yyref)) allocate(yyref(nres)) if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres) - if(.not.allocated(theta)) allocate(theta(nres)) - if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres)) - if(.not.allocated(theta)) allocate(theta(nres)) + if(.not.allocated(dc_norm)) then +! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2)) + allocate(dc_norm(3,0:2*nres+2)) + dc_norm(:,:)=0.d0 + endif + call int_from_cart(.true.,.false.) - call sc_loc_geom(.true.) -! wczesbiej bylo false + call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) enddo +! do i=1,2*nres +! vbld_inv(i)=0.d0 +! vbld(i)=0.d0 +! enddo + do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) enddo enddo - do i=2,nres-1 do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) @@ -2745,18 +3433,19 @@ write(iout,*) "nloctyp",nloctyp ! enddo ! enddo ! - allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm) - allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym) - allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym) + if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm) + if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym) + if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym) !----------------------------- 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) if (i.gt.1) then - if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then + if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then chain_length=lll-1 kkk=kkk+1 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) @@ -2789,6 +3478,9 @@ write(iout,*) "nloctyp",nloctyp ! 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) @@ -2825,7 +3517,7 @@ write(iout,*) "nloctyp",nloctyp do kkk=1,nperm write (iout,*) "nowa struktura", nperm do i=1,nres - write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),& + write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),& cref(2,i,kkk),& cref(3,i,kkk),cref(1,nres+i,kkk),& cref(2,nres+i,kkk),cref(3,nres+i,kkk) @@ -2850,29 +3542,10 @@ write(iout,*) "nloctyp",nloctyp enddo enddo ishift_pdb=ishift -!--------------------- -! el reallocate array - do i=1,2*nres+2 - do kkk=1,nperm - cref_alloc(1,i,kkk)=cref(1,i,kkk) - cref_alloc(2,i,kkk)=cref(2,i,kkk) - cref_alloc(3,i,kkk)=cref(3,i,kkk) - enddo - enddo -!el deallocate(cref) -!el allocate(cref(3,2*nres+2,nperm)) !(3,maxres2+2,maxperm) - - do i=1,2*nres+2 - do kkk=1,nperm - cref(1,i,kkk)=cref_alloc(1,i,kkk) - cref(2,i,kkk)=cref_alloc(2,i,kkk) - cref(3,i,kkk)=cref_alloc(3,i,kkk) - enddo - enddo -!--------------------- + return end subroutine readpdb -#ifndef WHAM_RUN +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) !----------------------------------------------------------------------------- ! readrtns_CSA.F !----------------------------------------------------------------------------- @@ -2922,7 +3595,7 @@ write(iout,*) "nloctyp",nloctyp character(len=640) :: controlcard real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,& - + integer i nglob_csa=0 eglob_csa=1d99 @@ -2961,6 +3634,65 @@ write(iout,*) "nloctyp",nloctyp 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 + protein=index(controlcard,"PROTEIN").gt.0 + ions=index(controlcard,"IONS").gt.0 + nucleic=index(controlcard,"NUCLEIC").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) dccart=(index(controlcard,'CART').gt.0) overlapsc=(index(controlcard,'OVERLAP').gt.0) @@ -3060,6 +3792,25 @@ write(iout,*) "nloctyp",nloctyp 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 !----------------------------------------------------------------------------- @@ -3182,7 +3933,7 @@ write(iout,*) "nloctyp",nloctyp ! 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.) @@ -3295,21 +4046,27 @@ write(iout,*) "nloctyp",nloctyp 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),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 @@ -3600,10 +4357,9 @@ write(iout,*) "nloctyp",nloctyp !----------------------------------------------------------------------------- subroutine openunits - use energy_data, only: usampl + use MD_data, only: usampl use csa_data use MPI_data -! use MD use control_data, only:out1file use control, only: getenv_loc ! implicit real*8 (a-h,o-z) @@ -3671,6 +4427,8 @@ write(iout,*) "nloctyp",nloctyp ! 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) @@ -3685,6 +4443,19 @@ write(iout,*) "nloctyp",nloctyp 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') @@ -3695,6 +4466,9 @@ write(iout,*) "nloctyp",nloctyp ! 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') @@ -3719,6 +4493,25 @@ write(iout,*) "nloctyp",nloctyp ! 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) + open (itube,file=tubename,status='old',action='read') + call getenv_loc('IONPAR',ionname) + open (iion,file=ionname,status='old',action='read') + ! print *,"Processor",myrank," opened file ISIDEP" ! print *,"Processor",myrank," opened parameter files" #elif (defined G77) @@ -3728,6 +4521,9 @@ write(iout,*) "nloctyp",nloctyp ! 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) @@ -3744,6 +4540,23 @@ write(iout,*) "nloctyp",nloctyp 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) + open (itube,file=tubename,status='old') + call getenv_loc('IONPAR',ionname) + open (iion,file=ionname,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& readonly) @@ -3752,6 +4565,8 @@ write(iout,*) "nloctyp",nloctyp ! 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) @@ -3774,11 +4589,50 @@ write(iout,*) "nloctyp",nloctyp 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('SIDEPAR_SCBASE',sidename_scbase) + open (isidep_scbase,file=sidename_scbase,status='old',action='read') + call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase) + open (isidep_pepbase,file=pepname_pepbase,status='old',action='read') + call getenv_loc('SCPAR_PHOSPH',pepname_scpho) + open (isidep_scpho,file=pepname_scpho,status='old',action='read') + call getenv_loc('PEPPAR_PHOSPH',pepname_peppho) + open (isidep_peppho,file=pepname_peppho,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') + call getenv_loc('IONPAR',ionname) + open (iion,file=ionname,status='old',action='read') + #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) 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 @@ -3939,8 +4793,12 @@ write(iout,*) "nloctyp",nloctyp thetname(:ilen(thetname)) write (iout,*) "Rotamer parameter file : ",& rotname(:ilen(rotname)) +!el---- +#ifndef CRYST_THETA write (iout,*) "Thetpdb parameter file : ",& thetname_pdb(:ilen(thetname_pdb)) +#endif +!el write (iout,*) "Threading database : ",& patname(:ilen(patname)) if (lentmp.ne.0) & @@ -3954,7 +4812,6 @@ write(iout,*) "nloctyp",nloctyp subroutine readrst use geometry_data, only: nres,dc - use energy_data, only: usampl,iset use MD_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' @@ -3967,10 +4824,15 @@ write(iout,*) "nloctyp",nloctyp open(irest2,file=rest2name,status='unknown') read(irest2,*) totT,EK,potE,totE,t_bath - do i=1,2*nres + totTafm=totT +! do i=1,2*nres +! AL 4/17/17: Now reading d_t(0,:) too + do i=0,2*nres read(irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo - do i=1,2*nres +! do i=1,2*nres +! AL 4/17/17: Now reading d_c(0,:) too + do i=0,2*nres read(irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo if(usampl) then