! Read molecular data.
!
use energy_data
- use geometry_data, only:nres,deg2rad,c,dc
+ use geometry_data, only:nres,deg2rad,c,dc,nres_molec
use control_data, only:iscode
use control, only:rescode,setup_var,init_int_table
use geometry, only:alloc_geo_arrays
!el integer :: rescode
!el real(kind=8) :: x(maxvar)
character(len=320) :: controlcard !,ucase
- integer,dimension(nres) :: itype_pdb !(maxres)
- integer :: i,j,i1,i2,it1,it2
+ integer,dimension(nres,5) :: itype_pdb !(maxres)
+ integer :: i,j,i1,i2,it1,it2,mnum
real(kind=8) :: scalscp
!el logical :: seq_comp
+ write(iout,*) "START MOLREAD"
+ call flush(iout)
+ do i=1,5
+ nres_molec(i)=0
+ print *,"nres_molec, initial",nres_molec(i)
+ enddo
+
call card_concat(controlcard,.true.)
call reada(controlcard,'SCAL14',scal14,0.4d0)
call reada(controlcard,'SCALSCP',scalscp,1.0d0)
call reada(controlcard,'TEMP0',temp0,300.0d0) !el
call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
r0_corr=cutoff_corr-delt_corr
- call readi(controlcard,"NRES",nres,0)
- allocate(sequence(nres+1))
+ call readi(controlcard,"NRES",nres_molec(1),0)
+ call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
+ call readi(controlcard,"NRES_CAT",nres_molec(5),0)
+ nres=0
+ do i=1,5
+ nres=nres_molec(i)+nres
+ enddo
+! allocate(sequence(nres+1))
!el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
call alloc_geo_arrays
call alloc_ener_arrays
!----------------------------
allocate(c(3,2*nres+2))
allocate(dc(3,0:2*nres+2))
- allocate(itype(nres+2))
+ allocate(itype(nres+2,5))
allocate(itel(nres+2))
+ if (.not. allocated(molnum)) allocate(molnum(nres+2))
!
! Zero out tableis.
do i=1,2*nres+2
enddo
enddo
do i=1,nres+2
- itype(i)=0
+ molnum(i)=0
+ do j=1,5
+ itype(i,j)=0
+ enddo
itel(i)=0
enddo
!--------------------------
!
+ allocate(sequence(nres+1))
iscode=index(controlcard,"ONE_LETTER")
if (nres.le.0) then
write (iout,*) "Error: no residues in molecule"
write (iout,*) "Error: too many residues",nres,maxres
endif
write(iout,*) 'nres=',nres
+! HERE F**** CHANGE
! Read sequence of the protein
if (iscode.gt.0) then
read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
endif
! Convert sequence to numeric code
- do i=1,nres
- itype(i)=rescode(i,sequence(i),iscode)
+ do i=1,nres_molec(1)
+ mnum=1
+ molnum(i)=1
+ itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
+ enddo
+ do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
+ mnum=2
+ molnum(i)=2
+ itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
+ enddo
+ do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
+ mnum=5
+ molnum(i)=5
+ itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
enddo
+
write (iout,*) "Numeric code:"
- write (iout,'(20i4)') (itype(i),i=1,nres)
+ write (iout,'(20i4)') (itype(i,molnum(i)),i=1,nres)
do i=1,nres-1
+ mnum=molnum(i)
#ifdef PROCOR
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
+ if (itype(i,mnum).eq.ntyp1_molec(mnum) .or. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
#else
- if (itype(i).eq.ntyp1) then
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
#endif
itel(i)=0
#ifdef PROCOR
- else if (iabs(itype(i+1)).ne.20) then
+ else if (iabs(itype(i+1,mnum)).ne.20) then
#else
- else if (iabs(itype(i)).ne.20) then
+ else if (iabs(itype(i,mnum)).ne.20) then
#endif
itel(i)=1
else
enddo
write (iout,*) "ITEL"
do i=1,nres-1
- write (iout,*) i,itype(i),itel(i)
+ mnum=molnum(i)
+ write (iout,*) i,itype(i,mnum),itel(i)
enddo
+ write(iout,*)
call read_bridge
if (with_dihed_constr) then
nnt=1
nct=nres
- if (itype(1).eq.ntyp1) nnt=2
- if (itype(nres).eq.ntyp1) nct=nct-1
+ if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
+ if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
write(iout,*) 'NNT=',NNT,' NCT=',NCT
call setup_var
call init_int_table
write (iout,'(20i4)') (iss(i),i=1,ns)
write (iout,'(/a/)') 'Pre-formed links are:'
do i=1,nss
+ mnum=molnum(i)
i1=ihpb(i)-nres
i2=jhpb(i)-nres
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,mnum)
+ it2=itype(i2,mnum)
write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
- restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',&
+ restyp(it1,molnum(i1)),'(',i1,') -- ',restyp(it2,molnum(i2)),'(',i2,')',&
dhpb(i),ebr,forcon(i)
enddo
endif
character(len=16) :: key
integer :: iparm
!el real(kind=8) :: ip,mp
- real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,&
+ real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
res1
integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n
- integer :: nlobi,iblock,maxinter,iscprol
+ integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
!
! Body
!
wvdwpp=ww(16)
wbond=ww(18)
wsccor=ww(19)
+ wcatcat=ww(44)
+ wcatprot=ww(43)
endif
!
weights(17)=wbond
weights(18)=0 !scal14 !
weights(21)=wsccor
+ weights(43)=wcatprot
+ weights(44)=wcatcat
+
! el--------
call card_concat(controlcard,.false.)
call reads(controlcard,"SCPPAR",scpname_t,scpname)
open (iscpp,file=scpname_t,status='old')
rewind(iscpp)
+ call getenv_loc('IONPAR',ionname)
+ open (iion,file=ionname,status='old')
+ rewind(iion)
write (iout,*) "Parameter set:",iparm
write (iout,*) "Energy-term weights:"
do i=1,n_eneW
endif
enddo
#else
- read (ibond,*) ijunk,vbldp0,akp,rjunk
+ read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
do i=1,ntyp
read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
j=1,nbondterm(i))
'inertia','Pstok'
write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
do i=1,ntyp
- write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
+ write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
vbldsc0(1,i),aksc(1,i),abond0(1,i)
do j=2,nbondterm(i)
write (iout,'(13x,3f10.5)') &
enddo
enddo
endif
+ if (.not. allocated(msc)) allocate(msc(ntyp1,5))
+ if (.not. allocated(restok)) allocate(restok(ntyp1,5))
+
+ 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
+
+
!----------------------------------------------------
allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
allocate(sig0(-ntyp:ntyp),sigc0(-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))
do i=1,ntyp
do j=1,ntyp
augm(i,j)=0.0D0
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)
+ write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp)
endif
! goto 50
!----------------------- LJK potential --------------------------------
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),&
+ write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
i=1,ntyp)
endif
! goto 50
read (isidep,*)(sigii(i),i=1,ntyp)
read (isidep,*)(chip(i),i=1,ntyp)
read (isidep,*)(alp(i),i=1,ntyp)
+ do i=1,ntyp
+ read (isidep,*)(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
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),&
+ write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
chip(i),alp(i),i=1,ntyp)
endif
! goto 50
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),&
+ write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
sigii(i),chip(i),alp(i),i=1,ntyp)
endif
case default
! 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(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
+ aa_aq(i,j)=0.0D0
+ bb_aq(i,j)=0.0D0
+ aa_lip(i,j)=0.0d0
+ bb_lip(i,j)=0.0d0
chi(i,j)=0.0D0
sigma(i,j)=0.0D0
r0(i,j)=0.0D0
do i=2,ntyp
do j=1,i-1
eps(i,j)=eps(j,i)
+ epslip(i,j)=epslip(j,i)
enddo
enddo
do i=1,ntyp
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)
if (ipot.gt.2) then
sigt1sq=sigma0(i)**2
sigt2sq=sigma0(j)**2
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,molnum(i)),restyp(j,molnum(i)),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
!-----------------------------------------------------------------------------
subroutine read_general_data(*)
- use control_data, only:indpdb,symetr
+ use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions
use energy_data, only:distchainmax
+ use geometry_data, only:boxxsize,boxysize,boxzsize
! implicit none
! include "DIMENSIONS"
! include "DIMENSIONS.ZSCOPT"
call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
call readi(controlcard,"RESCALE",rescale_modeW,1)
+ call reada(controlcard,'BOXX',boxxsize,100.0d0)
+ call reada(controlcard,'BOXY',boxysize,100.0d0)
+ call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+ ions=index(controlcard,"IONS").gt.0
+ 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
check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
call readi(controlcard,'SYM',symetr,1)
use control_data, only:pdbref
use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
nstart_sup,nstart_seq,nperm,nres0
- use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype
+ use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
use compare, only:seq_comp !,contact,elecont
use geometry, only:chainbuild,dist
use io_config, only:readpdb
character(len=4) :: sequence(nres)
!el integer rescode
!el real(kind=8) :: x(maxvar)
- integer :: itype_pdb(nres)
+ integer :: itype_pdb(nres,5)
!el logical seq_comp
- integer :: i,j,k,nres_pdb,iaux
+ integer :: i,j,k,nres_pdb,iaux,mnum
real(kind=8) :: ddsc !el,dist
integer :: kkk !,ilen
!el external ilen
return 1
34 continue
do i=1,nres
- itype_pdb(i)=itype(i)
+ mnum=molnum(i)
+ itype_pdb(i,mnum)=itype(i,mnum)
enddo
call readpdb
do i=1,nres
- iaux=itype_pdb(i)
- itype_pdb(i)=itype(i)
- itype(i)=iaux
+ iaux=itype_pdb(i,mnum)
+ itype_pdb(i,mnum)=itype(i,mnum)
+ itype(i,mnum)=iaux
enddo
close (ipdbin)
do kkk=1,nperm
nstart_seq=nnt
if (nsup.le.(nct-nnt+1)) then
do i=0,nct-nnt+1-nsup
- if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),&
+ if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
nsup)) then
do j=nnt+nsup-1,nnt,-1
do k=1,3
return 1
else
do i=0,nsup-(nct-nnt+1)
- if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),&
+ if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
nct-nnt+1)) &
then
nstart_sup=nstart_sup+i
enddo
enddo
do i=1,nres
+ mnum=molnum(i)
do j=1,3
dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
enddo
- if (itype(i).ne.10) then
+ if (itype(i,mnum).ne.10) then
ddsc = dist(i,nres+i)
do j=1,3
dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
use geometry_data, only:nres,c
- use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype
+ use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
'D','E','F','G','H','I','J'/),shape(chainid))
integer,dimension(nres) :: ica !(maxres)
real(kind=8) :: temp,efree,etot,entropy,rmsdev
- integer :: ii,i,j,iti,ires,iatom,ichain
+ integer :: ii,i,j,iti,ires,iatom,ichain,mnum
write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
ii,temp,rmsdev
write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
ichain=1
ires=0
do i=nnt,nct
- iti=itype(i)
+ mnum=molnum(i)
+ iti=itype(i,mnum)
if (iti.eq.ntyp1) then
ichain=ichain+1
ires=0
ires=ires+1
iatom=iatom+1
ica(i)=iatom
- write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
+ write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
ires,(c(j,i),j=1,3)
if (iti.ne.10) then
iatom=iatom+1
- write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
+ write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
ires,(c(j,nres+i),j=1,3)
endif
endif
enddo
write (ipdb,'(a)') 'TER'
do i=nnt,nct-1
- if (itype(i).eq.ntyp1) cycle
- if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
+ mnum=molnum(i)
+ if (itype(i,mnum).eq.ntyp1) cycle
+ if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
write (ipdb,30) ica(i),ica(i+1)
- else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
+ else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
write (ipdb,30) ica(i),ica(i+1),ica(i)+1
- else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
+ else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
write (ipdb,30) ica(i),ica(i)+1
endif
enddo
- if (itype(nct).ne.10) then
+ if (itype(nct,molnum(nct)).ne.10) then
write (ipdb,30) ica(nct),ica(nct)+1
endif
do i=1,nss
projects / unres4.git / blobdiff