! 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
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),&
!---------------------- 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)
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
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)
+ epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
+ alphi_scpho(j)
END DO
allocate(aa_scpho(ntyp_molec(1)))
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)
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
! 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(:)
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+1))
+ 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)
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
projects / unres4.git / blobdiff