X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_config.f90;h=3f4303ac0a17a5e96f77b87baf4ce8af642b1d10;hb=705644e0cbb7678faefd6fe1bc436159d38ad85d;hp=4c7f3a5971ac0b0276036b517d60c39779357f6e;hpb=e1aeda758c35ccdf5a5b88e825b6ca391c2dc1ae;p=unres4.git diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 4c7f3a5..3f4303a 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -735,11 +735,11 @@ 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,& @@ -782,9 +782,6 @@ allocate(nbondterm(ntyp)) !(ntyp) allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp) - 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)) @@ -792,6 +789,10 @@ 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) + read (ibond,*) vbldp0,akp,mp,ip,pstok do i=1,ntyp nbondterm(i)=1 @@ -804,6 +805,10 @@ endif enddo #else + 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),& @@ -854,6 +859,19 @@ 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) @@ -1412,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. @@ -1683,7 +1719,7 @@ #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) +! 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) @@ -1701,12 +1737,12 @@ !el-------------------- read (itorp_nucl,*,end=113,err=113) & (itortyp_nucl(i),i=1,ntyp_molec(2)) - print *,itortyp_nucl(:) +! 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) +! print *,nterm_nucl(i,j),nlor_nucl(i,j) v0ij=0.0d0 si=-1.0d0 do k=1,nterm_nucl(i,j) @@ -1866,6 +1902,7 @@ #endif if (lprint) then + l=3 write (iout,'(/a/)') 'Torsional constants:' do i=1,nsccortyp do j=1,nsccortyp @@ -1910,88 +1947,88 @@ do i=0,nloctyp-1 read (ifourier,*,end=115,err=115) - read (ifourier,*,end=115,err=115) (b(ii),ii=1,13) + read (ifourier,*,end=115,err=115) (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,'(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 @@ -2055,6 +2092,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)) @@ -2164,6 +2202,7 @@ end select continue close (isidep) + !----------------------------------------------------------------------- ! Calculate the "working" parameters of SC interactions. @@ -2272,6 +2311,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 @@ -2291,6 +2409,174 @@ 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 @@ -2344,6 +2630,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 ! @@ -2438,7 +2738,7 @@ ! Read the PDB file and convert the peptide geometry into virtual-chain ! geometry. use geometry_data - use energy_data, only: itype + use energy_data, only: itype,istype use control_data use compare_data use MPI_data @@ -2465,7 +2765,7 @@ character(len=80) :: card real(kind=8),dimension(3,20) :: sccor integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode, - integer :: isugar,molecprev + integer :: isugar,molecprev,firstion character*1 :: sugar real(kind=8) :: cou real(kind=8),dimension(3) :: ccc @@ -2475,6 +2775,7 @@ real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm) real(kind=8),dimension(:,:), allocatable :: c_temporary integer,dimension(:,:) , allocatable :: itype_temporary + integer,dimension(:),allocatable :: istype_temp efree_temp=0.0d0 ibeg=1 ishift1=0 @@ -2484,15 +2785,27 @@ ! 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) - + if(.not. allocated(istype)) allocate(istype(maxres)) do i=1,100000 read (ipdbin,'(a80)',end=10) card ! write (iout,'(a)') card @@ -2586,7 +2899,7 @@ ishift1=ishift1-1 !!!!! ! write (iout,*) "New chain started",ires,ishift,ishift1,"!" ires=ires-ishift+ishift1 - print *,ires,ishift,ishift1 +! print *,ires,ishift,ishift1 ires_old=ires ibeg=0 else @@ -2601,11 +2914,21 @@ 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 molecule=2 - itype(ires,molecule)=rescode(ires,res(2:4),0,molecule) + 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 @@ -2635,7 +2958,7 @@ .or. atom.eq."C4'" .or. atom.eq."O4'")) then read(card(31:54),'(3f8.3)') (ccc(j),j=1,3) !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3) - print *,ires,ishift,ishift1 +! print *,ires,ishift,ishift1 counter=counter+1 ! iii=iii+1 ! do j=1,3 @@ -2644,29 +2967,40 @@ 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(1:3).eq."C5'") then + else if (atom.eq."C5'") then read (card(19:19),'(a1)') sugar isugar=sugarcode(sugar,ires) if (ibeg.eq.1) then istype(1)=isugar else istype(ires)=isugar +! print *,ires,istype(ires) endif if (unres_pdb) then + 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. & @@ -2685,7 +3019,16 @@ 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 @@ -2697,7 +3040,9 @@ 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) + 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 @@ -2707,8 +3052,10 @@ ! 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(:) + 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 @@ -2725,7 +3072,7 @@ if (itype(i,k).eq.ntyp1_molec(k)) then if (itype(i+1,k).eq.ntyp1_molec(k)) then if (itype(i+2,k).eq.0) then - print *,"masz sieczke" +! print *,"masz sieczke" do j=1,5 if (itype(i+2,j).ne.0) then itype(i+1,k)=0 @@ -2740,47 +3087,47 @@ ! 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 +! 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 +! 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 +! endif !unres_pdb else !itype(i+1,1).eq.ntyp1 - if (unres_pdb) then +! 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 +! 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 +! 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 !unres_pdb endif !itype(i+1,1).eq.ntyp1 endif !itype.eq.ntyp1 @@ -2805,24 +3152,24 @@ if (molecule.eq.5) molecule=molecprev itype(nres,molecule)=ntyp1_molec(molecule) nres_molec(molecule)=nres_molec(molecule)+1 - if (unres_pdb) then +! 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)-1.9d0*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))/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(:) @@ -2946,7 +3293,10 @@ ! NOW LETS ROCK! SORTING allocate(c_temporary(3,2*nres)) allocate(itype_temporary(nres,5)) - allocate(molnum(nres)) + 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 @@ -2958,6 +3308,8 @@ 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 @@ -2971,6 +3323,7 @@ do k=1,5 do i=1,nres itype(i,k)=itype_temporary(i,k) + istype(i)=istype_temp(i) enddo enddo if (itype(1,1).eq.ntyp1) then @@ -3008,7 +3361,7 @@ enddo endif - print *,seqalingbegin,nres +! print *,seqalingbegin,nres if(.not.allocated(vbld)) then allocate(vbld(2*nres)) do i=1,2*nres @@ -3087,6 +3440,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) @@ -3579,7 +3933,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.) @@ -3692,21 +4046,27 @@ if(me.eq.king.or..not.out1file) & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",& eta - gamp=scal_fric*(pstok(i)+rwat)*eta - stdfp=dsqrt(2*Rb*t_bath/d_time) - allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1) +! allocate(gamp + do j=1,5 !types of molecules + gamp(j)=scal_fric*(pstok(j)+rwat)*eta + stdfp(j)=dsqrt(2*Rb*t_bath/d_time) + enddo + allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1) + do j=1,5 !types of molecules do i=1,ntyp - gamsc(i)=scal_fric*(restok(i,1)+rwat)*eta - stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) + gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta + stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time) enddo + enddo + if(me.eq.king.or..not.out1file)then write (iout,'(/2a/)') & "Radii of site types and friction coefficients and std's of",& " stochastic forces of fully exposed sites" - write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp) + write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1)) do i=1,ntyp write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),& - gamsc(i),stdfsc(i)*dsqrt(gamsc(i)) + gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1)) enddo endif else if (tbf) then @@ -4149,6 +4509,8 @@ 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" @@ -4193,6 +4555,8 @@ 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) @@ -4236,17 +4600,39 @@ 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