X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_config.f90;h=87ddf94beea67828febd369586e764f506f5e616;hb=10689ab7d813dfbdbb0c6e631d90234d78ea306a;hp=6186282609c4481372c539b431f5329aaf2d4f95;hpb=d62d8388175880fb6700daba86794663aeb3a78d;p=unres4.git diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 6186282..87ddf94 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -529,9 +529,9 @@ ! read secondary structure prediction from JPRED here! ! read(isecpred,'(A80)',err=100,end=100) line ! read(line,'(f10.3)',err=110) ftors - read(isecpred,'(f10.3)',err=110) ftors + read(isecpred,'(f10.3)',err=110) ftors(1) - write (iout,*) 'FTORS factor =',ftors + write (iout,*) 'FTORS factor =',ftors(1) ! initialize secstruc to any do i=1,nres secstruc(i) ='-' @@ -553,6 +553,7 @@ ii=0 do i=1,nres + ftors(i)=ftors(1) if ( secstruc(i) .eq. 'H') then ! Helix restraints for this residue ii=ii+1 @@ -735,7 +736,7 @@ 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 @@ -782,9 +783,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 +790,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 +806,15 @@ endif enddo #else + mp(:)=0.0d0 + ip(:)=0.0d0 + msc(:,:)=0.0d0 + isc(:,:)=0.0d0 + + 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),& @@ -1897,6 +1908,7 @@ #endif if (lprint) then + l=3 write (iout,'(/a/)') 'Torsional constants:' do i=1,nsccortyp do j=1,nsccortyp @@ -1941,88 +1953,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 @@ -2146,6 +2158,8 @@ !---------------------- GB or BP potential ----------------------------- case(3:4) ! 30 do i=1,ntyp +! print *,"I AM in SCELE",scelemode + if (scelemode.eq.0) then do i=1,ntyp read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp) enddo @@ -2173,6 +2187,81 @@ write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),& chip(i),alp(i),i=1,ntyp) endif + else +! print *,ntyp,"NTYP" + allocate(icharge(ntyp1)) +! print *,ntyp,icharge(i) + icharge(:)=0 + read (isidep,*) (icharge(i),i=1,ntyp) + print *,ntyp,icharge(i) +! if(.not.allocated(eps)) allocate(eps(-ntyp +!c write (2,*) "icharge",(icharge(i),i=1,ntyp) + allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp)) + allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp)) + allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp)) + allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp)) + allocate(epsintab(ntyp,ntyp)) + allocate(dtail(2,ntyp,ntyp)) + allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp)) + allocate(wqdip(2,ntyp,ntyp)) + allocate(wstate(4,ntyp,ntyp)) + allocate(dhead(2,2,ntyp,ntyp)) + allocate(nstate(ntyp,ntyp)) + if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1)) + if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp) + do i=1,ntyp + do j=1,i +! write (*,*) "Im in ALAB", i, " ", j + read(isidep,*) & + eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & + (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & + chis(i,j),chis(j,i), & + nstate(i,j),(wstate(k,i,j),k=1,4), & + dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& + dtail(1,i,j),dtail(2,i,j), & + epshead(i,j),sig0head(i,j), & + rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & + alphapol(i,j),alphapol(j,i), & + (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j) +! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) + END DO + END DO + DO i = 1, ntyp + DO j = i+1, ntyp + eps(i,j) = eps(j,i) + sigma(i,j) = sigma(j,i) + sigmap1(i,j)=sigmap1(j,i) + sigmap2(i,j)=sigmap2(j,i) + sigiso1(i,j)=sigiso1(j,i) + sigiso2(i,j)=sigiso2(j,i) +! print *,"ATU",sigma(j,i),sigma(i,j),i,j + nstate(i,j) = nstate(j,i) + dtail(1,i,j) = dtail(1,j,i) + dtail(2,i,j) = dtail(2,j,i) + DO k = 1, 4 + alphasur(k,i,j) = alphasur(k,j,i) + wstate(k,i,j) = wstate(k,j,i) + alphiso(k,i,j) = alphiso(k,j,i) + END DO + + dhead(2,1,i,j) = dhead(1,1,j,i) + dhead(2,2,i,j) = dhead(1,2,j,i) + dhead(1,1,i,j) = dhead(2,1,j,i) + dhead(1,2,i,j) = dhead(2,2,j,i) + + epshead(i,j) = epshead(j,i) + sig0head(i,j) = sig0head(j,i) + + DO k = 1, 2 + wqdip(k,i,j) = wqdip(k,j,i) + END DO + + wquad(i,j) = wquad(j,i) + epsintab(i,j) = epsintab(j,i) +! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j) + END DO + END DO + endif ! goto 50 !--------------------- GBV potential ----------------------------------- case(5) @@ -2201,9 +2290,11 @@ ! Calculate the "working" parameters of SC interactions. !el from module energy - COMMON.INTERACT------- - allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) + allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) + if (.not.allocated(chi)) allocate(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) + if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1)) + allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1) allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),& dcavtub(ntyp1)) allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),& @@ -2212,9 +2303,11 @@ bb_aq(:,:)=0.0D0 aa_lip(:,:)=0.0D0 bb_lip(:,:)=0.0D0 + if (scelemode.eq.0) then chi(:,:)=0.0D0 sigma(:,:)=0.0D0 r0(:,:)=0.0D0 + endif acavtub(:)=0.0d0 bcavtub(:)=0.0d0 ccavtub(:)=0.0d0 @@ -2230,6 +2323,7 @@ epslip(i,j)=epslip(j,i) enddo enddo + if (scelemode.eq.0) then do i=1,ntyp do j=i,ntyp sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2) @@ -2238,6 +2332,7 @@ rs0(j,i)=rs0(i,j) enddo enddo + endif if (lprint) write (iout,'(/a/10x,7a/72(1h-))') & 'Working parameters of the SC interactions:',& ' a ',' b ',' augm ',' sigma ',' r0 ',& @@ -2247,6 +2342,7 @@ epsij=eps(i,j) if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then rrij=sigma(i,j) +! print *,"SIGMA ZLA?",sigma(i,j) else rrij=rr0(i)+rr0(j) endif @@ -2257,6 +2353,7 @@ sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) aa_aq(i,j)=epsij*rrij*rrij + print *,"ADASKO",epsij,rrij,expon bb_aq(i,j)=-sigeps*epsij*rrij aa_aq(j,i)=aa_aq(i,j) bb_aq(j,i)=bb_aq(i,j) @@ -2268,7 +2365,7 @@ 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 + if ((ipot.gt.2).and. (scelemode.eq.0)) then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 sigii1=sigii(i) @@ -2435,7 +2532,7 @@ read(isidep_scbase,*) & dhead_scbasei(i,j), & dhead_scbasej(i,j), & - rborn_scbasei(i,j),rborn_scbasej(j,i) + 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) @@ -2500,6 +2597,77 @@ 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 @@ -2603,6 +2771,31 @@ write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,& ' v3ss:',v3ss endif + if (shield_mode.gt.0) then + pi=4.0D0*datan(1.0D0) +!C VSolvSphere the volume of solving sphere + 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*(4.50d0)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/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) +! if (scelemode.eq.0) then + short_r_sidechain(i)=sigma(i,i)/sqrt(2.0) + if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2 +! else +! short_r_sidechain(i)=sigma(i,i) +! endif + write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& + sigma0(i) + enddo + buff_shield=1.0d0 + endif + return 111 write (iout,*) "Error reading bending energy parameters." goto 999 @@ -2665,7 +2858,7 @@ use control_data use compare_data use MPI_data - use control, only: rescode,sugarcode +! use control, only: rescode,sugarcode ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.LOCAL' @@ -2708,18 +2901,30 @@ ! 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 + write (iout,'(a)') card if (card(:5).eq.'HELIX') then nhfrag=nhfrag+1 lsecondary=.true. @@ -2829,7 +3034,7 @@ ! 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) @@ -2899,11 +3104,19 @@ ! 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. & @@ -2944,7 +3157,8 @@ molecule=5 nres_molec(molecule)=nres_molec(molecule)+1 print *,"HERE",nres_molec(molecule) - itype(ires,molecule)=rescode(ires,res(2:4),0,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 @@ -2954,9 +3168,10 @@ ! Calculate dummy residue coordinates inside the "chain" of a multichain ! system nres=ires - if ((ires_old.ne.ires).and.(molecule.ne.5)) & + if (((ires_old.ne.ires).and.(molecule.ne.5)) & + ) & nres_molec(molecule)=nres_molec(molecule)-2 -! print *,'I have', nres_molec(:) + print *,'I have',nres, nres_molec(:) do k=1,4 ! ions are without dummy if (nres_molec(k).eq.0) cycle @@ -2988,47 +3203,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 +! 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 +! 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 @@ -3053,24 +3268,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(:) @@ -3172,10 +3387,10 @@ if (lprn) then write (iout,'(/a)') & "Cartesian coordinates of the reference structure" - write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & + write (iout,'(a,16x,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)') & + write (iout,'(5(a3,1x),i5,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 @@ -3186,7 +3401,7 @@ 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)') & + write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') & ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),& (c(j,nres+ires),j=1,3) enddo @@ -3194,7 +3409,9 @@ ! 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 @@ -3225,36 +3442,36 @@ istype(i)=istype_temp(i) enddo enddo - if (itype(1,1).eq.ntyp1) then - nsup=nsup-1 - nstart_sup=2 - if (unres_pdb) then +! 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 +! 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))') & + write (iout,'(a,16x,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)') & + write (iout,'(5(a3,1x),i5,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 @@ -3342,7 +3559,7 @@ write (iout,*) "symetr", symetr do i=1,nres lll=lll+1 -!c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) +! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) if (i.gt.1) then if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then chain_length=lll-1 @@ -3371,7 +3588,7 @@ ! write (iout,*) "spraw lancuchy",chain_length,symetr ! do i=1,4 ! do kkk=1,chain_length -! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3) +! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3) ! enddo ! enddo ! enddiagnostic @@ -3394,7 +3611,7 @@ cou=0 do kkk=1,symetr icha=tabperm(i,kkk) -! write (iout,*) i,icha + write (iout,*) i,icha do lll=1,chain_length cou=cou+1 if (cou.le.nres) then @@ -3402,7 +3619,7 @@ kupa=mod(lll,chain_length) iprzes=(kkk-1)*chain_length+lll if (kupa.eq.0) kupa=chain_length -! write (iout,*) "kupa", kupa + write (iout,*) "kupa", kupa cref(j,iprzes,i)=chain_rep(j,kupa,icha) cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha) enddo @@ -3421,11 +3638,11 @@ cref(3,i,kkk),cref(1,nres+i,kkk),& cref(2,nres+i,kkk),cref(3,nres+i,kkk) enddo - 100 format (//' alpha-carbon coordinates ',& + 100 format (//' alpha-carbon coordinates ',& ' centroid coordinates'/ & ' ', 6X,'X',11X,'Y',11X,'Z', & 10X,'X',11X,'Y',11X,'Z') - 110 format (a,'(',i3,')',6f12.5) + 110 format (a,'(',i5,')',6f12.5) enddo !c enddiag @@ -3555,6 +3772,13 @@ call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) call readi(controlcard,'TUBEMOD',tubemode,0) + print *,"SCELE",scelemode + call readi(controlcard,"SCELEMODE",scelemode,0) + print *,"SCELE",scelemode + +! elemode = 0 is orignal UNRES electrostatics +! elemode = 1 is "Momo" potentials in progress +! elemode = 2 is in development EVALD write (iout,*) TUBEmode,"TUBEMODE" if (TUBEmode.gt.0) then call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) @@ -3692,22 +3916,22 @@ write (iout,'(2a)') diagmeth(kdiag),& ' routine used to diagonalize matrices.' if (shield_mode.gt.0) then - pi=3.141592d0 + pi=4.0D0*datan(1.0D0) !C VSolvSphere the volume of solving sphere -!C print *,pi,"pi" + 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 + VSolvSphere=4.0/3.0*pi*(4.50d0)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/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 +! 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 @@ -3867,6 +4091,10 @@ large = index(controlcard,"LARGE").gt.0 print_compon = index(controlcard,"PRINT_COMPON").gt.0 rattle = index(controlcard,"RATTLE").gt.0 + preminim=(index(controlcard,'PREMINIM').gt.0) + write (iout,*) "PREMINIM ",preminim + dccart=(index(controlcard,'CART').gt.0) + if (preminim) call read_minim ! if performing umbrella sampling, fragments constrained are read from the fragment file nset=0 if(usampl) then @@ -4503,6 +4731,11 @@ 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')