! 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) ='-'
ii=0
do i=1,nres
+ ftors(i)=ftors(1)
if ( secstruc(i) .eq. 'H') then
! Helix restraints for this residue
ii=ii+1
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
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))
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
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),&
#endif
if (lprint) then
+ l=3
write (iout,'(/a/)') 'Torsional constants:'
do i=1,nsccortyp
do j=1,nsccortyp
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
!---------------------- 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
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)
! 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),&
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
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)
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 ',&
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
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)
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)
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)
do i=1,ntyp_molec(1)
do j=1,ntyp_molec(2)-1
epsij=eps_scbase(i,j)
- rrij=sigma(i,j)
+ rrij=sigma_scbase(i,j)
! r0(i,j)=rrij
! r0(j,i)=rrij
rrij=rrij**expon
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
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
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'
! 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.
! 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)
! 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. &
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
! 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
! 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
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(:)
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
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
! 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
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
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
! 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
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
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
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
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)
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
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
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')
projects / unres4.git / blobdiff